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Abstract 

A  rotor  with  internal  damping  is  a  complex  dynamic  system.  An  accurate  model  is 
required  to  design  and  test  controllers  for  supercritical  operation.  This  thesis  developed 
and  validated  a  model  for  this  purpose. 

An  existing  rotordynamic  Finite  Element  Method  (FEM)  model  and  magnetic 
bearing  simulation,  developed  by  Draper  Laboratory,  were  first  compared  to  actual  rotor 
test  data.  Correlation  of  predicted  and  actual  parameters  such  as  critical  speeds,  rigid 
body  modes,  and  rotor  displacements  was  used  to  validate  the  Draper  model  and 
simulation.  Once  this  correlation  was  established,  both  the  model  and  simulation  were 
modified  to  include  the  destabilizing  effects  of  internal  damping.  Using  these  improved 
predictive  tools,  several  Proportional  Integral  and  Derivative  (PID)  controllers  were 
designed  and  implemented  in  an  effort  to  stabilize  a  rotor  system  with  internal  damping. 
The  PID  controllers  were  effective  in  stabilizing  the  rotor  systems  at  sub-critical  speeds. 
However,  the  model  developed  during  this  thesis  showed  that  these  controllers  were 
unable  to  counteract  the  destabilizing  effects  of  internal  damping  during  supercritical 
operation.  This  improved  rotordynamic  model  and  magnetic  bearing  simulation  is  now 
available  to  test  more  complex  controller  designs  in  the  supercritical  regime. 
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Nomenclature 


A  -  Cross  sectional  area  of  dashpot  damper 

[A]  -  State  space  dynamic  matrix 

[B]  -  State  space  actuator  matrix 

Ci,C2  -  Amplitudes  of  single  degree  of  freedom  solution 

c  -  Viscous  damping  constant 

cc  -  Critical  damping  constant 

ce  -  External  viscous  damping  constant 

c;  -  Internal  viscous  damping  constant 

[Cd]  -  External  damping  matrix,  [Cd]  =  2[£c d] 

[C]  -  State  space  sensor  matrix 

[Ce]  -  External  damping  matrix 

D  -  Energy  dissipated  per  cycle 

E  -  Modulus  of  Elasticity 

E'  -  Loss  Modulus 

{f}  -  Fixed  frame  external  forces 

F  -  Forcing  function  amplitude 

F,x  -  Force,  displacement  (structural) 

[G]  -  Gyroscopic  matrix 

[G]  -  Normalized  gyroscopic  matrix 

A 

[G]  -  Normalized  gyroscopic,  external  damping  matrix  (internal  damping  only) 

G  -  Shear  Modulus 

h  -  Height  of  dashpot  damper 

i  -  Complex  quantity  V-I 

I  -  Moment  of  Inertia 

Id  -  Diametral  Moment  of  Inertia 

Ip  -  Polar  Moment  of  Inertia 

Im  [  ]  -  Imaginary  part  of  [  ] 

[7]  -  Identity  matrix,  diagonal 

k  -  Stiffness 

ky  -  Elements  in  stiffness  matrix  [K] 

[K]  -  Total  stiffness  matrix  (includes  bending  and  shear  deformations) 

[Kbend]  -  Bending  stiffness  matrix 

[KciJ  -  Circulation  stiffness  matrix 

[KShear]  -  Shear  stiffness  matrix 

Kd  -  Derivative  Constant  for  PID  Controller 

Ki  -  Integral  Constant  for  PID  Controller 

Kp  -  Proportional  Constant  for  PID  Controller 

[L]  -  fix  m  matrix  that  inflates  the  force  matrix,  {f } ,  from  an  m  x  1  to  an  h  x  1  matrix 

l  -  Element  length 

m  -  Mass 

m  -  Number  of  forces  applied  to  the  rotor 
M  -  Moment 


[M]  -  Mass  matrix 

[Mtransl  -  Translational  mass  matrix 

[Mfot]  -  Rotational  mass  matrix 

M',V'  -  Moment  and  shears  on  the  right  side  of  the  stations  in  the  Myklestad-Prohl 
method 

h  -  Number  of  modes  in  modal  matrix 
Pn  -  Elements  of  the  modal  matrix  (undamped  eigenvectors) 
pi,p2  -  Exponential  part  of  single  degree  of  freedom  solution 
{q}  -  Fixed  frame  physical  coordinates 

r  -  Radial  location  of  element  on  shaft,  internal  viscous  damping 
r0  -  Radius  of  shaft,  viscous  internal  damping 

R  -  Radius  of  whirl,  internal  viscous  damping 

Re[  ]  -  Real  part  of  [  ] 

Sj  -  Complex  eigenvalue  solution 

t  -  Time 

T  -  Kinetic  energy 

u  -  Velocity  in  dashpot  damper 

V  -  Potential  energy 

v  -  Velocity  in  dashpot  damper 

V  -  Centerline  displacement  in  y  direction 

V*,W”  -  Second  derivatives  of  displacements 
Vshear  -  Contribution  to  displacement  due  to  shear 
Vbend  -  Contribution  to  displacement  due  to  bending 
<W>  -  Time  average  power  dissipated 

W  -  Centerline  displacement  in  x  direction 

x  -  Displacement 

x  -  Amplitude  of  phasor 

x,x,x  -  Displacement,  velocity,  acceleration 

y  -  Distance  from  center  to  edge  of  dashpot  damper 

y0  -  Shaft  deflection,  internal  viscous  damping 

Zm  -  Mechanical  impedance 
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{a}  -  State  space  coordinates 

(3  -  Rotation  about  y-axis  (Jeffcott  Rotor)  or  Rotation  around  x  axis  (Nelson) 

8  -  Logarithmic  decrement  or  Rotation  about  x-axis,  Jeffcott  Rotor 

e  -  Strain 

e0  -  Strain  at  edge  of  shaft,  internal  viscous  damping 

i  -  Logarithmic  decrement  for  the  internal  hysteretic  damping  of  the  shaft  divided 

by  n 

<j)  -  Angular  shaft  displacement  around  x-axis  (hysteretic) 

r|  -  Loss  factor 

[T|]  -  Internal  damping  matrix 

Tjeff  -  Effective  loss  factor 

TlH-struct  ‘  Structural  hysteretic  loss  factor 

T|h  -  Internal  hysteretic  damping  coefficient 

T)v  -  Internal  viscous  damping  coefficient 

(p  -  Phase  angle 

[<J>]  -  Modal  matrix 

m  -  Normalized  modal  matrix 

[d>]  -  Shear  deformation 

r  -  Rotation  around  y-axis 

k  -  cross  sectional  shape  factor  for  shear  deformation  of  a  shaft 

A  -  Translational  shape  functions  for  Timoshenko  rotor  element 

A,  -  Whirl  orbit  grow/decay 

p  -  Viscosity 

Pi  -  Viscosity,  internal  damping 

pe  -  Viscosity,  external  damping 

0  -  Angle  between  shaft  whirl  and  shaft  rotation,  internal  viscous  damping 

0  -  Angular  shaft  displacement  around  y-axis  (hysteretic) 

0  -  Rotational  shape  functions  for  Timoshenko  rotor  element 

p  -  Density 

G  -  Stress 

G  -  Stress  (complex  quantity) 

d  -  Stress  under  sinusoidal  motion 

x  -  Shear  stress 

to  -  Whirl  speed 

COd  -  Frequency  of  damped  vibration 

con  -  Natural  frequency 

[CD2]  -  Diagonal  matrix,  squared  eigenvalues 

Q  -  Rotor  spin  speed 

Qcr  -  Critical  Speed 

£2th  -  Instability  Threshold 

{£}  -  Modal  coordinates 

£  -  Damping  ratio 
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Chapter  1  Introduction 


1.1  Background 

The  rotors  in  many  commercial  and  military  rotating  machines  operate  above  the 
first,  second  and  even  third  critical  speeds,  [1].  Substantial  weight  and  performance 
advantages  are  afforded  to  those  machines  that  can  run  at  these  high  speeds,  [2].  But  as 
the  speed  increases,  eventually  the  rotor  will  become  unstable. 

A  primary  destabilizing  mechanism  for  these  high-speed  rotors  is  the  presence  of 
internal  damping.  This  type  of  damping,  due  to  strain  of  the  shaft  material  or  to  micro¬ 
movements  at  shrink  fits  and/or  couplings,  can  under  certain  circumstances,  feed  energy 
into  the  shaft  vibration  and  thereby  induce  instability,  [2]. 

Recently,  magnetic  bearings  have  emerged  as  a  possible  solution  to  the 
destabilizing  effects  of  internal  damping.  The  basic  advantage  of  magnetic  bearings  over 
standard  bearings  is  that  the  stiffness  and  damping  of  the  magnetic  bearings  can  be  varied 
as  a  function  of  spin  speed  or  other  conditions  of  the  rotating  system.  This  advantage  has 
already  been  exploited  by  the  development  of  closed  loop  controllers  that  increase  the 
damping  of  the  magnetic  bearings  as  the  rotor  passes  through  a  critical  speed,  and  then 
decrease  the  damping  after  critical  speed  passage. 

Furthermore,  as  the  spin  speed  approaches  an  instability  threshold,  it  should  be 
possible  to  develop  a  controller  to  increase  the  external  damping  of  the  magnetic  bearings 
(or  make  other  modifications)  to  counteract  the  destabilizing  effects  of  the  rotor’s  internal 
damping.  This  increased  external  damping  (or  other  factors)  should  be  able  to  delay  the 
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onset  of  instability  and  allow  the  rotor  to  be  run  at  even  higher  speeds.  This  thesis  will 
clarify  the  concepts  of  internal  damping;  validate  an  existing  FEM  rotor  model  and 
magnetic  bearing  simulation  by  comparing  predictions  to  actual  rotor  test  data; 
incorporate  internal  damping  into  the  rotor  model  and  simulation;  and  it  will  test  a  basic 
PID  controller  to  determine  if  it  is  sufficient  to  counteract  the  effect  of  internal  damping. 

1.2  Thesis  Motivation 

In  line  with  its  work  on  Robust  Control  Schemes  for  Magnetic  Bearings,  The 
Charles  Stark  Draper  Laboratory  has  developed  a  rotordynamic  finite  element  model 
(Antkowiak)  that  can  successfully  predict  critical  speeds.  In  addition,  the  program  also 
produces  a  state  space  model  of  the  rotordynamic  system  for  use  in  a  magnetic  bearing 
simulation  (Scholten).  Although  the  model  and  simulation  were  highly  detailed,  they  did 
not  incorporate  any  destabilizing  effects.  By  adding  internal  damping  to  this  FEM 
model,  the  instability  threshold  of  the  system  can  be  predicted,  yielding  a  more  accurate 
state  space  model  for  supercritical  operation.  In  addition,  this  model  could  now  be 
placed  in  a  magnetic  bearing  simulation  so  that  different  controllers  could  be  designed 
and  tested  to  achieve  improved  control  and  possibly  increase  the  instability  threshold, 
thus  allowing  the  rotor  to  run  at  higher  speeds. 

In  addition  to  the  control-based  goals  mentioned  above,  the  implementation  of 
state-of-the-art  technology  into  current  Draper  Laboratory  projects  is  an  equally 
important  goal  of  this  thesis.  The  most  appropriate  project  for  the  magnetic  bearing  and 
control  technology  is  the  Draper  Laboratory  Flywheel  Energy  Storage  effort.  This 
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Internal  Research  and  Development  project  is  devoted  to  demonstrating  the  first  flywheel 
energy  storage  system  capable  for  use  in  space.  This  thesis  is  aimed  at  contributing  a 
high  fidelity  model  and  improved  simulation  to  assist  in  the  testing  of  various  controller 
schemes  that  could  allow  the  flywheel  to  operate  at  higher  speeds;  thus  allowing  the 
storage  of  more  energy.  The  ability  to  predict  unstable  rotordynamic  behavior  and  give  a 
certain,  limited  amount  of  control  to  counteract  these  effects  can  also  benefit  many  other 
ongoing  and  future  projects  in  the  Laboratory. 

1.3  Thesis  Outline 

This  chapter  briefly  discusses  some  of  the  challenges  in  the  field  of  rotordynamics 
and  introduces  the  possibility  of  using  magnetic  bearings  to  compensate  for  the 
destabilizing  effects  of  internal  damping.  The  motivations  for  the  thesis  are  then 
described  and  the  goals  are  presented.  A  brief  overview  of  subsequent  chapters  follows. 

Chapter  Two  discusses  the  concepts  of  external  and  internal  damping.  First, 
external  damping  is  described,  including  the  external  viscous  damping  model.  Internal 
damping  is  then  defined  and  further  divided  into  two  distinct  types:  internal  viscous  and 
internal  hysteretic  damping.  The  models  for  each  case  are  presented  in  detail. 

Chapter  Three  takes  the  concept  of  internal  damping  and  applies  it  to  rotating 
systems.  Again  the  distinction  between  internal  viscous  and  internal  hysteretic  damping 
is  made.  An  emphasis  is  placed  on  the  role  of  internal  damping  in  determining  the 
instability  thresholds.  Next,  the  specific  effect  of  external  damping  on  the  instability 
threshold  is  addressed. 


1-3 


Chapter  Four  introduces  the  finite  element  method  as  a  means  to  analyze  complex 
rotordynamic  systems.  The  internal  damping  models  are  defined,  and  the  element 
matrices  and  the  equations  of  motion  are  referenced  and  summarized.  The  model  is  then 
validated  using  a  standard  test  case  from  the  references. 

Chapter  Five  discusses  the  development  of  a  detailed  state  space  model 
(Antkowiak)  and  a  high  fidelity  simulation  (Scholten)  that  includes  internal  damping. 
First  the  method  that  the  FEM  program  used  to  produce  the  state  space  model  is  outlined, 
then  the  full  simulation  configuration  is  presented.  Finally,  some  details  are  given  on  the 
specific  PID  controller  that  was  used  in  the  simulation. 

Chapter  Six  presents  the  original  FEM  modeling  and  simulation  predictions,  for  a 
rotor  without  internal  damping,  and  compares  them  to  actual  rotor  test  data.  First,  the 
rotor  FEM  model  is  described  for  an  actual  rotor  located  at  Draper  Laboratory.  Next,  the 
model  predictions  for  the  rigid  body  modes  and  critical  speeds  are  given.  Using  the  state 
space  model  generated  by  the  FEM  program,  the  simulation  predictions  for  the  rigid  body 
modes,  critical  speeds,  and  rotor  displacements  are  presented.  Finally,  the  FEM  and 
simulation  predictions  are  compared  to  actual  rotor  test  results. 

Chapter  Seven  presents  the  modeling  and  simulation  predictions  for  a  rotor  with 
various  combinations  of  internal  damping.  First,  instability  threshold  predictions  from 
the  FEM  program  are  presented.  Then,  the  full  simulation  results  are  given,  and  the 
effectiveness  of  the  PID  controller  is  discussed. 

The  thesis  is  summarized  and  conclusions  are  drawn  in  Chapter  Eight. 
Recommendations  are  made  for  future  work  at  Draper  Laboratory  and  also  for  the  fields 
of  rotordynamics  and  controls  in  general,  using  magnetic  bearing  technology. 
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Chapter  2  Damping 


When  a  shaft  spins,  it  experiences  a  jump  rope  like  motion  called  whirl  or 
precession,  due  to  slight  shaft  imbalance  or  other  external  forces.  This  whirling  motion 
can  put  the  shaft  into  alternating  states  of  stress/strain.  Since  the  shaft  material  is  not 
perfectly  elastic,  work  is  performed  on  the  system  with  every  rotation.  This  additional 
energy  added  to  the  whirling  rotor  system  is  due  to  the  internal  damping  of  the  shaft.  As 
the  speed  and  work  done  on  the  shaft  increases,  the  whirl  orbit  of  the  shaft  begins  to 
grow.  At  lower  speeds,  the  bearings  of  the  rotating  system  are  able  to  dissipate  this 
energy/work  by  transmitting  it  to  the  ground/surroundings  (external  damping).  However, 
as  the  rotation  speed,  energy,  and  shaft  amplitude/orbit  increase,  the  bearings  will  not  be 
able  to  dissipate  all  of  the  energy  and  the  shaft  will  become  unstable. 

There  are  two  common  models  of  internal  damping:  viscous  and  hysteretic. 
Viscous  internal  damping  is  defined  as  a  velocity  dependent  quantity  that  is  represented 
as  a  dashpot  with  coefficient  T|y  This  type  of  internal  damping  tends  to  destabilize  the 
system  only  after  the  first  critical  speed.  Hysteretic  internal  damping  is  defined  as  a 
displacement  dependent  quantity  expressed  by  the  loss  factor  T|h.  This  type  of  damping 
can  be  destabilizing  at  any  speed  (except  for  zero  speed). 

2.1  External  Damping  (Viscous  Model) 

Prior  to  elaborating  on  internal  damping,  the  more  familiar  “external  damping”,  or 
damping  to  ground,  will  be  briefly  discussed.  If  a  system  with  only  mass  and  elasticity, 
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completely  isolated  from  its  surroundings,  is  disturbed  by  a  perturbation,  it  will  begin  to 
vibrate.  Since  the  system  is  not  connected  to  its  surroundings,  no  energy  can  be 
transferred  away  from  the  system  and  it  will  continue  to  vibrate  forever.  Since  this 
situation  never  occurs  in  nature,  an  energy  transfer  mechanism  must  exist  between  the 
system  and  its  surroundings,  [3].  In  many  systems,  this  transfer  mechanism  gradually 
converts  vibrational  energy  to  heat  or  sound.  This  is  known  as  damping  or  external 
damping,  since  the  energy  is  taken  from  the  vibrating  system  and  transferred  to  the 
“external”  surroundings,  [4]. 

The  most  widely  used  model  for  external  damping  is  the  linear  viscous  model. 
This  model  assumes  that  energy  is  removed  from  a  system  in  the  same  manner  as  energy 
is  removed  by  the  use  of  a  dashpot.  A  dashpot  is  a  piston-cylinder  type  device  filled  with 
viscous  oil  or  other  fluid.  As  the  piston  is  moved  down  into  the  fluid  (by  a  force,  with  an 
associated  energy)  the  fluid  passes  by  the  piston  in  a  thin  fluid  layer  (see  Figure  1). 
According  to  Newton’s  Laws  of  viscous  flow,  the  shear  stress  developed  in  the  fluid 
layer  between  the  piston  and  the  wall  of  the  cylinder  is  given  by 


du 


v 

h 


(1) 


The  shear  resisting  force  developed  along  the  piston  is 


-cv 


(2) 


where 


pA 

h 


is  defined  as  the  viscous  damping  constant  “c”. 


Therefore,  the  Newtonian 


viscous  damping  model  predicts  a  linear  relationship  between  the  damping  force  and  the 
velocity.  More  specifically,  the  damper  exerts  a  force  that  is  proportional  to  and  in  anti¬ 
phase  with  the  velocity,  [4]. 
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Surface  Area,  A 


Figure  1:  Newtonian  Viscous  Damping  Model 

In  a  single  degree  of  freedom  system,  the  viscous  external  damping  is  represented 
by  the  cx  term  in  the  standard  equation  of  motion  of 

mx  +  cx  +  kx  =  FcosQt  (3) 

By  assuming  the  solution  to  be  of  the  form  x  =  ept ,  the  homogeneous  equation  (with 
right-hand  side  equal  to  zero)  can  be  solved  to  give: 

x(t)  =  Clep'‘  +C2ePi‘ ,  (4) 


where 


and  where  the  natural  frequency  of  the  system  is  given  by 


The  critical  damping,  ec,  is  defined  as  the  value  of  the  damping  constant  for  which 
the  terms  under  the  radical  in  equation  (5)  become  zero.  Further,  the  damping  ratio  is 
defined  as 


and  the  frequency  of  damped  vibration  is  given  by 


(8) 


ad=o)J  1-C2 

Using  all  of  these  definitions,  the  homogenous  solution  or  transient  solution  to  equation 
(3)  takes  the  form  of: 

*  =  (9) 

If  the  damping  ratio  is  less  than  one,  underdamped  oscillation  occurs.  This  is 
characterized  by  oscillatory  motion  and  by  the  amplitude  of  vibration  decreasing 
exponentially  with  time.  If  the  damping  ratio  is  equal  to  one,  the  system  is  critically 
damped.  In  this  case,  the  response  of  the  system  will  be  aperiodic  and  will  quickly 
diminish  to  zero  without  oscillation.  Finally,  if  the  damping  ratio  is  greater  than  one,  the 
system  is  overdamped  and  the  amplitude  will  also  diminish  to  zero  without  oscillation 
(although  not  quite  as  fast  as  the  critically  damped  case),  [4], 

2.2  The  Linear  Viscoelastic  Model 

Viscoelastic  materials  are  materials  which  exhibit  both  elastic  and  viscous 
behavior,  [13].  A  true  representative  of  a  viscoelastic  material  is  a  polymer.  In  this  type 
of  material,  the  atoms  are  joined  strongly  together,  although  the  long  chains  can  be 
branched  or  physically  entangled  (providing  either  a  weak  or  strong  link).  While  under  a 
cyclic  load,  the  elastic  nature  of  the  polymer  arises  from  the  stretching  of  the 
intermolecular  bonds  or  entangled  chains.  In  addition,  the  breaking  and  reforming  of  the 
intermolecular  bonds  contributes  to  the  viscous  nature,  or  damping,  of  the  material,  [5]. 

When  dealing  with  forced  (steady  state)  sinusoidal  motions,  it  is  convenient  to 
represent  the  stress  and  strain  by  rotating  vectors  in  the  complex  plane.  Typically,  these 
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vectors  are  called  phasors.  For  the  case  of  a  material  under  a  sinusoidal  motion  of 
frequency  Q,  the  stress  can  be  given  as 

G  =  Oi (10) 

where  a  represents  a  complex  quantity.  Similarly,  the  strain  of  the  material  under  the 
motion  is  given  by 

e=eeiat  (11) 

The  physical  stress  and  strain  can  be  found  by  taking  the  real  part  of  the  respective 
phasors: 

o(t)  =  Re[a]  and  e(t)  =  Re[f  ] .  (12) 

To  show  the  relationship  between  the  stress  and  the  strain  of  a  material  under 
sinusoidal  load,  the  two  phasors  are  drawn  on  the  complex  plane  (Figure  2).  From 
observation,  it  is  has  been  found  that  the  strain  lags  the  stress  by  the  phase  angle  (p,  [6]. 


Figure  2:  Stress  and  Strain  Phasor  Diagram 

Plotting  this  relationship  on  the  real  stress-strain  axes  produces  the  classic 
hysteresis  ellipse  (Figure  3).  The  dashed  line  represents  a  system  with  no  losses,  or  the 
constitutive  relationship  by  Hooke  of  o  =  Ee .  As  the  losses  of  the  system  increase,  the 
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size  of  the  ellipse  also  begins  to  increase.  Based  on  the  definitions  of  work  and  energy, 
the  enclosed  area  of  the  ellipse  is  the  energy  dissipated  per  unit  volume  per  cycle  of 
oscillation,  [6, 7]. 


Figure  3:  Stress-Strain,  Energy  Dissipation  Ellipse  for  Viscoelastic  Material 

The  behaviors  in  Figures  2  and  3  can  be  summarized  in  the  equation: 

(13) 

where  E  represents  a  complex  modulus  of  the  form: 

E_  =  E  +  iE'.  (14) 

The  imaginary  part  E'  represents  the  loss  modulus  due  to  the  material  and  is  given  by 
E'  =  Etaxup .  (15) 

where  E  is  the  traditional  Modulus  of  Elasticity. 

Note  that  this  definition  of  the  total  complex  modulus  and  the  loss  modulus  are  consistent 
with  Figures  2  and  3. 

Substituting  the  complex  modulus  into  Hooke’s  constitutive  relationship  and 
factoring  out  E  gives: 


<7  — 


E' 

E{\  +  i—)e. 
E 


(16) 
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Now,  defining  the  loss  factor  r|  as 


E' 

fl  =  ~E  =tm<p’ 

gives  the  constitutive  law  of  a  one-dimensional,  linear  viscoelastic  material: 

o  =  E{\  +  iij)e. 

This  is  represented  by  the  model  in  Figure  4,  [1 1, 13, 17]. 


(17) 


(18) 


Figure  4:  Linear  Model  of  Viscoelastic  Material 

It  is  important  to  note  that  this  linear  model  is  the  basis  for  both  types  of  internal 
damping.  The  difference  between  viscous  internal  damping  and  hysteretic  internal 
damping  lies  in  how  the  loss  factor  term  is  defined.  Although  it  is  well  known  that  linear 
models  of  all  types  have  limitations  (the  specific  limitations  of  linear  internal  damping 
were  demonstrated  by  Graham,  [8]),  the  overriding  benefits  of  a  linear  stability 
evaluation  have  been  recognized  as  a  meaningful  indication  of  the  effects  of  internal 
damping  on  a  rotor  system,  [2]. 
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2.3  The  Internal  Viscous  Damping  Model 


One  of  the  simplest  mathematical  models  for  describing  internal  viscous  damping 
is  the  Voigt  Model.  It  takes  the  linear  viscoelastic  model  (shown  in  Figure  4)  and 
designates  the  damping  mechanism  as  an  idealized  dashpot  with  a  viscosity  of  |ij.  Figure 
5  shows  this  standard  model  for  viscous  internal  damping,  [7]. 


Figure  5:  Linear  Internal  Viscous  Damping  Model 


The  constitutive  equation  for  this  model  is  given  by: 

de 

o  =  Ee  +  v-, 


(19) 


or 


qeia>  =  Eee™  +  ju, 


int 


d(£e*“) 
dt 


(20) 


If  the  strain  derivative  is  evaluated  and  the  common  exponential  term  is  eliminated,  the 
equation  reduces  to: 

g_  =  (E  +  i£ljUi)£.  (21) 

Now  comparing  the  term  to  the  r\  term  in  equation  (18)  gives  the  effective  loss  factor 
of 
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ET]eff  =  /t,Q  or  I y  = 


(22) 


where  the  T|v  represents  the  internal  viscous  damping  coefficient  defined  as 

>7v=4-  <23) 

E 

This  result  indicates  that  the  effective  loss  factor  for  internal  viscous  damping  increases 
linearly  with  frequency,  [11, 13, 17]. 

The  internal  damping  coefficient  T|v  is  related  to  the  internal  damping  dashpot 
constant,  Ci,  if  one  inserts  equation  (23)  into  the  following  equation  derived  by  Gunter 
[9]: 


where  |lj  is  the  internal  viscosity  and  /  is  the  length  of  the  material.  Note  equation  (24)  is 
for  a  beam  (shaft)  only. 

To  determine  the  energy  dissipated  per  unit  cycle,  the  dashpot  in  the  viscous 
internal  damping  model  is  treated  exactly  the  same  as  an  external  viscous  dashpot.  This 
yields  the  familiar  equation  of  motion: 

mx+c^  +  kx- FcosQt  (25) 

where  again  q  is  the  internal  damping  constant.  Similar  to  an  external  damper,  the 
internal  damper  exerts  a  force  that  is  proportional  to  and  in  anti-phase  with  the  velocity. 
The  mechanical  impedance  of  equation  (25)  is  given  by: 

Zm=c,+i(mQ—^)  (26) 

and  the  time  average  power  dissipated  by  an  oscillator  is  given  by  the  equation: 
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<W>=^R  e[ZJ 


(27) 


where  the  Jc2  term  is  the  amplitude  of  the  velocity  phasor  or  Qx ,  [5,  6].  To  calculate  the 
energy  dissipated  per  cycle  the  following  equation  is  used: 

D=<w>^,[5]  (28) 

The  energy  dissipated  per  cycle  for  the  viscous  damping  case  is  then  calculated  to  be: 

D  =  7iQx2ci  (29) 

Like  equation  (21),  this  also  confirms  that  the  internal  viscous  damping  model  is  rate 
dependent;  in  particular,  it  is  linear  with  frequency,  [4,  6]. 

Another  consequence  of  this  rate  dependence  can  be  seen  from  Figure  3.  Since 
the  angle  (p  in  Figure  3  is  related  to  the  loss  factor,  from  equation  (17),  both  depend  on 
the  frequency.  This  means  that  as  the  frequency  changes,  the  size  of  the  ellipse  also 
changes. 

2.4  The  Internal  Hysteretic  Damping  Model 

The  first  researches  to  identify  internal  hysteretic  damping  in  rotating  systems 
were  Newkirk  and  Kimball  in  1924  [8,  10].  They  conducted  rotordynamic  tests  on  a 
series  of  different  materials  over  a  frequency  range  of  0.03  to  50  Hz  at  low  stress 
amplitudes.  Their  experiment  involved  connecting  a  weight  to  the  “overhang”  end  of  the 
shaft  by  a  bearing  as  shown  in  Figure  6.  This  allowed  the  shaft  to  rotate  freely  with  a 
constant  downward  force.  Newkirk  and  Kimball  found  that  as  the  shaft  was  rotated,  a 
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small  deflection  of  the  weight  to  one  side  occurred  no  matter  what  rotation  speed  was 
used. 

After  an  analysis,  they  concluded  that  the  angular  deflection  cp,  shown  in  Figure  6, 
was  cause  by  a  force  associated  with  the  internal  hysteretic  damping  of  the  material.  By 
measuring  this  angle  cp,  and  using  it  to  calculate  the  amount  of  work  done  by  the  system 
to  overcome  the  internal  damping/friction,  the  log  decrement  associated  with  the  material 
was  determined.  Since  these  initial  tests,  much  experimental  evidence  has  been  gathered 
to  confirm  that  internal  hysteretic  damping  in  a  material  is  rate-independent,  [7,  10].  In 
addition,  it  has  been  found  that  this  damping  force  is  in  anti-phase  with  the  velocity  but 
proportional  to  the  displacement,  [6]. 


Figure  6:  Kimball/Newkirk  Testing  Apparatus 

To  find  the  material  loss  factor  from  the  measured  angle  (p,  Kimball  and  Newkirk 
used  the  relationships  between  the  stress  and  strain  derived  in  Section  2.2.  A  summary  is 
shown  in  Figure  7. 
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Figure  7:  Stress  and  Strain  Phasor  Diagram  (Modified) 

Their  analysis  proceeds  as  follows:  starting  from  the  basic  sine  relationship  of 
E' 

sin  (p  = 


(30) 


Factoring  out  E  from  the  denominator  and  using  equation  (17)  gives: 

r 

E'  E  _  Vh 


sin  <p  =  - 


eM  + 


\E  j 


1+ 


'e')2 

l¥J 


•y/r+77/T 


(31) 


Hence  if  tp  is  measured  by  the  apparatus  of  Figure  6,  T|h  can  be  solved  from  the  quadratic 
relationship  given  by  equation  (31).  Again,  for  these  types  of  materials,  E  and  T|  are 
essentially  constant.  In  contrast  to  the  viscous  model,  having  a  constant  loss  factor 
dictates  that  the  energy  dissipation  ellipse  in  Figure  3  does  not  change  size  with  changing 
frequencies.  It  is  also  important  to  point  out  that  the  relationship  between  the  material 
loss  factor  T|h  and  the  structural  loss  factor  (riH-stmct )  depends  on  the  geometry  and 
loading  of  the  structure.  [6] 

To  calculate  the  axial  stress/strain  constitutive  equation  for  the  hysteretic  model 
shown  in  Figure  8,  the  components  of  the  complex  storage  modulus  Fare  used. 
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Figure  8:  Linear  Internal  Hysteretic  Damping  Model 

The  component  along  the  real  axis  in  Figure  7  (associated  with  Young’s  Modulus)  is 
found  using  the  familiar  in-phase  stress  strain  of: 

Re[o  ]  =  Ecos  (pe  (32) 

The  component  in  the  imaginary  direction  of  the  complex  storage  modulus  is  found  from 
the  quadrature  stress-strain  law  which  gives: 

Im[<7]  =  _Esin^>—  (33) 


Together,  they  give 


&  =  E\ 


ecos  (p-\ — sin  (p 


(34) 


Recalling  equation  (31)  and  recognizing  that 
1 


cos  (p  = 


V1  +  77tf2 


gives  the  constitutive  relationship 


<7  =  E 


e  +.  £Vh 


V1+77w2  W1  +  7?//2 


which  agrees  with  references  [2,  8]. 


(35) 


(36) 
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Writing  equation  (36)  in  complex  notation  gives: 


(Je‘n'  =  E\ 


£e 


iQt 


* 


+  -{ee'a,)X  Vh 


+  71h 


2  dt 


Q 


(37) 


Evaluating  the  derivative  and  canceling  out  the  appropriate  terms  yields: 


cr  =  E\ 


1  •  1 H 

r  +  l-  'H 


(38) 


or 


ar  =  E- 


From  Figure  7,  it  can  be  shown  that 


(39) 


(40) 


If  this  substitution  is  made,  the  resulting  equation 
0_  =  E(l  +  iT}H)£ 


(41) 


is  in  agreement  with  equation  (18)  and  the  previous  results  given  in  section  2.2. 

It  is  important  to  note  that  unlike  equation  (21),  the  frequency,  Q,  no  longer 
remains  in  equation  (41)  when  the  material  has  internal  hysteretic  damping,  validating  the 
experimental  results  obtained  by  Kimball  and  Newkirk. 

To  calculate  the  energy  dissipated  per  unit  cycle,  one  must  transform  the  variables 
from  local  material  properties  like  E  (Young’s  Modulus)  to  global  structural  properties 
like  k  (stiffness),  to  give: 

F  =  k(\  +  irjH_stmc,)x.  (42) 


2-14 


In  this  equation  the  force  and  displacements  are  written  as  phasors  while  the  structural 
stiffness  is  given  by  k,  [6]. 

Now,  if  the  lumped  mass  is  included  in  the  system,  equation  (42)  becomes: 

F  =  mx  +  k(l  +  irjH)x  .  (43) 

The  mechanical  impedance  of  this  equation  is  found  to  be: 

Z„=^  +  «mfl-£)  (44) 

Using  the  time  average  power  dissipation  equation,  equation  (27),  the  energy  dissipated 
per  unit  cycle  for  the  hysteretic  model  is  calculated  to  be: 

D  =  7d2kJ]H  (45) 

This  also  confirms  that  the  internal  hysteretic  damping  model  is  independent  of 
frequency,  [6]. 

Although,  the  hysteretic  damping  model  is  typically  associated  with  internal 
damping,  occasionally,  some  have  used  this  model  for  external  damping.  This  external 
hysteretic  damping  model  has  the  same  frequency  independence  and  energy  dissipation 
characteristics  as  the  internal  model.  But  unlike  internal  damping,  the  external  hysteretic 
damping  model  is  always  connected  physically  to  the  surrounding  and  always  dissipates 
energy  from  the  system,  [6]. 
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2.5  Internal  Viscous  Damping  vs.  Internal  Hysteretic  Damping 


Once  again,  the  difference  between  internal  viscous  damping  and  internal 
hysteretic  damping  lies  in  how  the  loss  factor  term  is  defined.  For  internal  viscous 
damping,  the  r|v  coefficient  has  been  identified  with  the  damping  constant  in  equation 
(24),  and  was  also  shown  to  increase  linearly  with  frequency,  [6, 7,  8].  This  damping 
force  associated  with  the  c(. x  term  is  proportional  to  and  in  anti-phase  with  the  velocity. 
For  the  internal  hysteretic  damping  coefficient,  r|H,  experimental  and  analytical  methods 
have  been  used  to  confirm  its  rate-independence,  [10].  This  damping  force  is  in  anti¬ 
phase  with  the  velocity  but  proportional  to  the  displacement. 

To  get  a  clear  graphical  explanation  of  the  differences,  recall  that  the  enclosed 
area  of  the  hysteresis  ellipse  in  Figure  3  is  the  energy  dissipated  per  cycle.  For  viscous 
internal  damping,  the  angle  between  the  stress  and  strain  increases  as  the  frequency 
increases,  resulting  in  a  larger  and  larger  ellipse.  In  contrast,  since  the  hysteretic  internal 
damping  coefficient  and  the  angle  remain  constant,  the  ellipse  does  not  change  size  or 
orientation  with  changing  frequency. 

To  summarize  the  differences  and  to  show  how  the  models  compare  to  some 
examples  of  real  materials.  Figure  9  and  Table  1  are  provided. 
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D 


Figure  9:  Variation  of  Energy  Dissipated  per  Cycle  with  Frequency,  [5] 


Table  1:  Forced  Excitation  of  a  Single  Degree  of  Freedom  System  with  Viscous  or 
_  Hysteretic  Damping,  [5]  _ 


Viscous  Damping 

Hysteretic  Damping 

Differential  Equation 

mx  +  cx  +  kx  =  FcosQt 

mx + k(  1  +  ir\H  )x  =  Re[Fe'n'  ] 

Particular  Solution 

F 

F 

■yj(k-mo)2)2  +(cco)2 

yJ(k-mco2)2  +(kTjH)2 

Energy  dissipated/cycle 

D  =  nQx2c 

D  =  nx2kT)H 

Resonant  Frequency 

Decreases  with  increasing 
value  of  c 

Independent  of  value  of  r)H 

Static  Displacement  at  x  =  0 

F 

k 

Depends  on 

Resonant  Amplitude 

Depends  on  all  equation 
parameters 

Independent  of  mass 

From  the  figure  and  reference  [5],  it  is  clear  that  most  systems  do  not  appear  to  have  a 
substantial  amount  of  viscous  damping.  Although  most  materials  do  exhibit  some 
characteristics  of  hystereric  internal  damping,  many  experts  admit  that  the  subject  matter 
is  not  well  understood,  [5,  8,  1 1].  Nevertheless,  combinations  and  slight  modifications 
of  the  two  models  appear  to  be  useful  for  many  applications. 
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Chapter  3  -  Internal  Damping  in  Rotors 


3.1  Rotordynamic  Basics 

The  simplest  rotor  model  consists  of  a  rigid  circular  disc  centered  on  a  massless 
shaft  supported  by  rigid  bearings.  The  rigid  disc  is  only  allowed  to  move  in  the  x-y 
plane.  This  model,  the  Jeffcott  or  Laval  rotor,  is  shown  in  Figure  1 1 .  The  free  vibration 
of  this  rotor  is  very  similar  to  the  two  degree  of  freedom  system  in  Figure  12  where  the 
spring  and  dashpot  pairs  are  perpendicular  to  one  another.  The  equations  of  motion  for 
any  rotation  speed  are  given  by: 

mx+cx+k c  =  0 

(46) 

my +  cy  +  ky  =  0 

The  solution  to  these  two  equations  describes  the  path  of  the  shaft  center  and  does  not 
depend  on  the  rotation  speed  of  the  shaft.  In  general,  the  natural  motion  of  the  shaft  takes 
the  form  of  an  ellipse  or  in  special  cases  a  circle  or  a  straight  line,  [12]. 


Figure  11:  Jeffcott  Rotor 
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Figure  12:  Jeffcott  Rotor,  Model  for  Natural  Vibration 


If  the  circular  disc  is  not  centered  on  the  shaft,  then  gyroscopic  forces,  the 
moments  of  inertia  of  the  disc,  and  the  shaft  spin  speed  all  must  be  taken  into  account, 
[12].  A  Jeffcott  rotor  with  a  non-central  disc  has  equations  of  motion  given  by: 


mx+k22x-k  23/6  =  0 
my  +  k22  +  yk23S  =  0 
Id8  +  IpClft+k23y  +  k33S  =  0 
ij-ipns-k23x+k33fi= 0 


(47a,b,c,d) 


where  the  displacements  are  given  by  x  and  y  and  the  rotations  are  given  by  8  and  (3  as  in 


Figure  13.  The  diametral  and  polar  moments  of  inertia  are  given  by  Id  and  Ip  respectively 
while  elements  in  the  stiffness  matrix  are  given  by  ky.  The  moments  associated  with  the 
polar  moment  of  inertia  and  the  rotation  speed  in  equations  (47c,d)  are  gyroscopic 


moments. 


3-2 


Figure  13:  Jeffcott  Rotor,  Non-Centered  Disc 

By  writing  the  equations  in  (47)  in  matrix  notation,  and  assuming  an  exponential 
solution,  four  pairs  of  complex  conjugate  eigenvalues  result  in  the  form  Sj  =  ±i<Bj . 

These  four  natural  frequencies  result  from  the  gyroscopic  terms  in  equations  (47)  and  the 
fact  that  the  system  has  four  degrees  of  freedom.  A  Campbell  Diagram,  such  as  Figure 
14,  shows  how  the  eigenvalues  change  as  the  shaft  spin  speed  increases,  [12]. 


Figure  14:  Campbell  Diagram  for  Jeffcott  Rotor 
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Based  on  the  eigenvectors  or  mode  shapes  of  the  solution,  the  positive  roots  were 
labeled  as  the  forward  whirl  frequencies,  while  the  negative  roots  were  labeled  as  the 
backward  whirl  frequencies.  These  definitions  also  correspond  to  physical  observation  as 
the  positive  roots  were  demonstrated  by  forward  whirl  (precession  or  whirl  in  the  same 
direction  of  the  rotation  of  the  shaft)  while  the  negative  roots  resulted  in  backward  whirl 
(precession  in  the  opposite  direction  of  the  rotation  of  the  shaft),  [12]. 

On  Figure  14,  it  is  important  to  point  out  two  places  where  the  rotor  spin  speed 
matches  that  of  the  forward  natural  whirling  frequency  of  the  rotor.  These  two  locations 
are  called  critical  speeds.  A  critical  speed,  is  characterized  by  a  large  local  amplitude 
vibration,  similar  to  resonance,  [9]. 

An  important  difference  between  the  model  with  a  non-central  disc  (or  any 
shaft/rotor  with  gyroscopic  forces)  and  the  centered  disc  model  is  that  the  natural  motion 
of  the  non-central  disc  is  circular.  In  addition,  it  should  be  noted  that  the  Jeffcott  rotor 
with  the  centered  disc  only  has  a  single  natural  frequency  while  the  models  with 
gyroscopic  forces  have  four  which  depend  on  the  shaft  rotation  speed,  [12]. 

The  simple  model  of  a  rotor  with  a  non-centered  disc  can  be  further  extended  to 
include  many  other  factors.  Flexible  bearings  with  external  damping,  an  unbalance  in  the 
disc,  and  all  forms  of  internal  damping  in  the  shaft  can  also  be  included  to  improve  the 
model.  The  discussion  by  Kramer  in  [12]  is  a  good  reference  that  shows  the  development 
of  these  improved  models. 
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3.2  Interna]  Viscous  Damping 


In  1933,  Smith,  [6],  evaluated  the  effects  of  internal  viscous  damping  on  a 
rotating  system.  He  proved  that  by  including  internal  viscous  damping  in  a  rotating 
system  without  external  damping,  the  system  would  become  unstable  at  the  first  critical 
speed.  This  speed  was  called  the  instability  threshold.  He  also  showed  that  internal 
viscous  damping  actually  has  stabilizing  effects  if  the  system  is  kept  at  sub-critical 
speeds,  and  destabilizing  effects  above  the  critical  speed. 

A  few  years  later,  the  study  by  Smith  was  confirmed  by  a  number  of  authors. 
These  included,  but  were  not  limited  to,  Kramer,  [12],  Ehrich,  [13],  Gunter,  [9],  and 
Dimarogonas,  [14].  These  authors  clarified  Smith’s  claim  and  showed  that  the 
introduction  of  internal  viscous  damping  caused  unstable,  sub-synchronous  rotor  motion 
above  the  lowest  critical  speed,  never  below.  Although  the  methods  each  author  used  to 
calculate  the  actual  instability  threshold  differed,  all  of  the  general  conclusions  were  the 
same.  In  addition,  they  also  noted  that  if  external  damping  is  added  to  the  system,  the 
stability  threshold  can  be  made  greater  than  the  lowest  critical  speed,  [9]. 

A  direct  and  simple  method  to  achieve  the  same  results  was  introduced  by  Ehrich, 
[13].  He  first  began  with  the  constitutive  relationship  given  by  equation  (19): 


Using  the  linearized  beam  theory  or  plane  sections  remain  plane  (where  the  strain  is 


proportional  to  the  distance  from  the  neutral  axis)  yields: 
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'  Bearing  Center  Line 

Figure  15:  Rotating,  Whirling  Beam  Cross  Section,  Ehrich,  [13] 

In  this  figure,  the  shaft  is  whirling  about  the  bearing  centerline  with  a  speed  of  co 
and  spinning  about  the  center  of  the  shaft  with  speed  of  Q.  The  y-axis  is  an  extension  of 
the  whirl  radius  R  and  also  rotates  about  the  bearing  centerline  with  a  speed  (D.  Note  that 
the  z-axis  remains  perpendicular  to  the  y-axis  and  the  shaft  rotates  about  the  x-axis  with 
respect  to  the  y-z  plane.  Looking  at  the  angles  defined  by  the  figure, 

6  =Qt- at  or  6  =Q.-co  (49) 

where  6  is  the  rotation  rate  of  the  shaft  with  respect  to  x,  y,  and  z.  If  the  motion  is 
synchronous,  i.e.  £2  =  a) ,  then  6=0  and  the  shaft  does  not  rotate  with  respect  to  the  x,  y, 
and  z  axes.  Assuming  sub-synchronous  motion  and  combining  equations  (19),  (48),  and 
(49),  gives 


E  cos  6  -  //,  (£2  -  co)  sin  6] . 


(50) 


If  this  stress  is  multiplied  by  the  distance  from  the  neutral  axis  of  the  shaft  (z-axis)  and 


integrated  over  the  cross  sectional  area,  the  resulting  moments  become: 
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(51) 


M  =  [  crxrcos  OdA  =  —EI 
Ja  ro 

My=j(JxrsmddA  =  -^iiiI(n-(o)  (52) 

A  ro 

The  moment  in  the  z  direction  is  due  to  the  bending  of  the  shaft  in  the  x-y  plane 
(bearings  remain  fixed,  shaft  experiences  jump  rope  like  motion,  or  whirl).  Since  the 
shaft  is  not  bent  in  the  x-z  plane,  the  moment  in  the  y  direction  must  be  balanced  by 
additional  forces.  To  achieve  equilibrium,  a  viscous  internal  damping  force  must  be 
tangent  to  the  shaft  whirl  orbit  (see  Figure  16).  The  particular  direction  of  the  force 
depends  on  the  sign  of  the  (Q-co)  term.  For  sub-synchronous  motion  it  is  positive  (i.e.  the 
spin  speed  is  greater  than  the  whirl  speed).  In  this  case  the  force  is  in  the  direction  of  the 
whirling  shaft  and  feeds  energy  into  the  system,  [6]. 


if  £2  >  (0, 

^internal  viscous  damping  “  UNSTABLE 


Figure  16:  Rotating,  Whirling  Beam  Cross  Section 
Showing  Direction  of  Internal  Damping  Force,  No  External  Damping 


Recall  the  claim  that  any  amount  of  external  damping  can  improve  the  stability  of 
the  system.  This  fits  into  the  above  argument  due  to  the  fact  that  these  external  damping 
forces  are  always  in  the  opposite  direction  of  the  shaft  whirl,  thereby  removing  energy 
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from  the  system.  The  stability  lies  in  the  comparison  of  the  two  forces  (see  Figure  17).  If 
the  external  damping  force  dominates,  the  shaft  is  stable,  but  if  the  internal  damping  force 
dominates,  the  shaft  is  unstable.  The  instability  threshold  occurs  where  they  exactly 
cancel. 


(O,  F internal  viscous  damping 


Figure  17 :  Rotating,  Whirling  Beam  Cross  Section 
Showing  Directions  of  Internal  and  External  Damping  Forces 


Ehrich  continues  in  the  derivation  to  define  values  of  both  the  internal  and 
external  forces  and  obtains: 


Finiema, damping  =  //,  ~  G>)1 =  C(-  (H ~  CO)  (53) 

Vernal  damping  =  MeWo  =  Wo®  (54) 


where  y0  is  the  maximum  shaft  deflection. 

At  the  instability  threshold,  the  two  forces  are  set  equal  to  one  another  giving: 

ci(Q.-co)  =  ceco  (55) 
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At  this  point,  the  shaft  will  be  whirling  at  its  critical  speed,  (Qcr).  This  has  been 
confirmed  by  observations  made  by  Kimball,  [10],  and  by  a  separate  analysis  by  Ehrich, 
[13].  Making  theta  =  Qcr  substitution  and  rearranging  the  terms  in  a  non-dimensional 
form  leads  to: 


Q 


th  _ 


a 


f  c  ' 

1+-*- 

V  Ci  J 


(56) 


which  agrees  with  Smith’s  original  equation,  [8]. 

To  graphically  represent  the  effects  of  viscous  internal  damping,  Figure  18  shows 
a  typical  Campbell  Diagram  with  whirl  speed  ©  and  rotor  spin  speed  Q.  The  first  critical 
speed,  QCr,  and  the  instability  threshold  speed,  Qth,  are  shown  on  the  spin  speed  axis.  In 
this  case,  since  the  rotor  does  not  become  unstable  at  the  first  critical  speed,  a  certain 
amount  of  external  damping  exists.  Eventually,  as  the  rotation  speed,  energy,  and  shaft 
amplitude/orbit  increase,  the  bearings  will  not  be  able  to  dissipate  all  of  the  energy  and 
the  shaft  will  become  unstable. 


Rotor  Spin  Speed,  Q. 


Figure  18:  Campbell  Diagram  for  Rotor  with  Internal 
and  External  Viscous  Damping 
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3.3  Internal  Hysteretic  Damping 


As  noted  previously  in  Chapter  2,  the  first  researchers  to  study  internal  hysteretic 
damping  in  rotating  systems  were  Newkirk  and  Kimball  in  1924,  [10].  While  performing 
their  tests  and  doing  analyses,  they  concluded  that  internal  hysteretic  damping  might  be  a 
possible  cause  for  the  shaft  whirl  and  instability.  Years  later,  Dimentberg,  [15],  and 
Lund,  [11],  took  advantage  of  the  work  done  by  Kimball  and  incorporated  internal 
hysteretic  damping  into  the  rotordynamic  equations  of  motion. 

Lund  began  his  derivation  with  the  internal  hysteretic  damping  model  (defined  in 
Chapter  2).  The  constitutive  relationship  between  the  bending  strain  and  bending  stress 
for  this  model  was  originally  expressed  by  Lund  as: 


—  =  —M.  cos 7 — —M  sin  7  = - /  (m  -eM  ) 

dz  El  z  El  }  EI-J l  +  ¥ 

v  e  (57a,b) 

—  =  —  M7sinr+—  Afvcosr  = - /  •  (m  +iM  ) 

dz  El  z  El  >  Eli 1+i7 

where  0  and  <|)  are  angular  shaft  displacement  defined  in  Figure  19.  The  logarithmic 

decrement  for  the  internal  hysteretic  damping  of  the  shaft  was  related  to  the  quantity  i 

by 

(58) 


„  6 
£=  — 
n 


and 


S  = 


2nX 
a j 


(59) 
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where  the  X  and  the  to  are  the  real  and  imaginary  parts  of  the  eigenvalue  solution  to  the 
equations  of  motion.  Note  that  the  constitutive  equations  already  suggest  an  iterative 
solution. 


Figure  19:  Rotor  Configuration  for  Displacements,  Lund 

Note:  Later  Nelson,  [2],  refined  the  constitutive  relationship  into  a  more  familiar  form 
relating  stress  and  strain  directly,  as  described  in  Chapter  2: 


&  =  E\ 


+  W*+77h2  _ 


(36) 


Next,  Lund  modified  the  conventional  FEM  beam  equations  to  include  internal 
hysteretic  damping  and  then  converted  them  into  a  form  suitable  for  the  Myklestad-Prohl 
numerical  method.  This  method  involved  representing  the  rotor  by  a  series  of  lumped 
masses  located  at  stations  connected  by  uniform  shaft  sections,  (somewhat  similar  to  the 
modem  finite  element  method).  This  technique  produced  the  following  equations  used  to 
solve  for  the  displacement  at  the  next  station,  n+1,  based  on  the  displacement  of  the 
previous  one: 
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*»+l  =*n+lA  + 


yn+ 1  =  +  lnt  + 


1 


mj- I 

i 


+i2 


(£/)> 


+  £2 


t-(M;,+a/L)+ 


_ (i»), 

6  (*GA), 


(®), 


(60a, b) 


V'  -gp' ) 

6  (*GA)„  |V  ,B  J 


Mx,n+l  =M'xn+lnV:n 

My,n+1=M'yn+lnV;n 


V  =V' 

r  x,n+l  r  xn 

V  =v' 

'  )  ,n+1  vyn 


(61a,b) 


(62,a,b,c,d) 


where  the  moments  and  shears  are  represented  by  M  and  V,  and  the  shaft  material  and 
physical  properties  as  K,  G,  E,  A  and  I  (see  Figure  20). 


Station  n 


Station  n+1 


Figure  20:  Rotor  Configuration  for  Bending  and  Shear  Forces,  Lund 

In  addition  to  the  displacements  of  the  rotor,  the  Myklestad-Prohl  method  was 
also  used  to  solve  for  the  eigenvalues  of  the  equations  of  motion.  This  solution  took  the 
form  of  S  =  A  ±  ico ,  where  the  real  part  indicated  the  growth/decay  of  the  whirl  orbit  and 
the  imaginary  part  represented  the  whirl  speed  noted  previously.  Lund  then  further 
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reduced  this  result  into  a  more  manageable  form  by  the  use  of  the  log  decrement  given  in 
equation  (59): 

S  =  ~— -  (63) 

O) 

Lund  added  the  negative  sign  to  indicate  that  positive  values  of  the  logarithmic 
decrement  resulted  in  stable  systems  while  negative  values  resulted  in  unstable  systems. 

Using  this  log  decrement  notation,  the  case  for  an  undamped  (no  external 
damping)  rotor  system  with  internal  hystereric  damping  was  solved.  Lund  found  that  the 
log  decrement  for  all  forward  modes  was  always  negative,  which  indicated  instability.  In 
addition,  he  found  that  the  log  decrement  for  all  backward  modes  was  always  positive, 
which  indicted  stability.  In  a  separate  study  mentioned  by  Lund,  [1 1],  he  added  external 
damping  through  the  bearings  to  the  rotating  system,  and  found  that  the  forward  mode 
stability  could  be  achieved  for  all  speeds.  The  results  and  conclusions  obtained  by  Lund 
agreed  with  the  earlier  analysis  conducted  by  Dimentberg,  [15]. 

3.4  The  Effect  of  External  Damping  on  the  Instability  Threshold 

To  assist  in  the  study  of  the  effect  of  external  damping  on  the  instability 
threshold,  the  Draper  Laboratory  FEM  Rotor  program  was  employed  (see  Chapter  4  for 
discussion  and  verification  of  this  program).  For  this  short  numerical  study,  viscous 
external  damping  was  applied  separately  to  two  rotor  systems:  one  with  internal  viscous 
damping  and  one  with  internal  hysteretic  damping. 
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The  first  case  that  was  run  on  the  program,  verified  the  results  obtained  by  Smith 
and  Ehrich,  where  the  instability  threshold  for  a  system  with  internal  and  external  viscous 


damping  was  given  by: 


Q 


Ih  _ 


Q_ 


1-A 

V  C‘  J 


(56) 


Figure  21,  shows  the  non-dimensional  ratios  of  the  instability  threshold  and  the 
critical  speed  versus  the  ratio  of  the  external  viscous  damping  to  internal  viscous 
damping.  The  data  from  the  finite  element  model  solution  match  closely  with  equation 
(56).  As  more  finite  elements  were  added  to  the  program,  the  results  continued  to 
approach  the  exact  solution. 


Figure  21:  Stability  Thresholds  for  Viscous  Internal  and 
Viscous  External  Damping 


The  next  configuration  involved  adding  viscous  external  damping  to  a  rotor 
system  with  only  hysteretic  internal  damping.  Since  the  ratio  of  viscous  external 
damping  and  hysteretic  internal  damping  cannot  be  put  in  a  non-dimensional  form,  a  plot 
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different  from  Figure  21  was  required  to  show  the  relationship.  Figure  22  shows  the 
values  of  viscous  external  damping  (given  by  Jie/E)  required  to  maintain  stability  at  all 
speeds  for  given  values  of  internal  hysteretic  damping.  For  zero  internal  hysteretic 
damping,  no  viscous  external  damping  is  required  to  maintain  stability.  For  increasing 
values  of  internal  hysteretic  damping,  the  required  viscous  damping  to  maintain  stability 
at  all  speeds  increases  linearly.  If  the  external  viscous  damping  is  less  than  that  given  by 
Figure  22,  all  of  the  forward  modes  will  be  unstable  in  a  range  of  spin  speeds  beginning 
at  Q  =  0+. 


Figure  22:  Amount  of  Viscous  External  Damping  Required  for  Stability  with 
Internal  Hysteretic  Damping 

Another  way  to  demostrate  this  relationship  is  shown  in  Figure  23.  This  plots  the 
log  decrement  of  the  eigenvalue  solution  (as  defined  by  Lund,  [11])  for  fixed  internal 
hysteretic  damping  as  a  function  rotor  spin  speed  Q.  It  also  shows  a  family  of  curves, 
which  represent  constant  viscous  external  damping.  Although  it  seems  unlikely  that  the 
rotating  system  will  ever  be  spun  at  the  high  speeds  shown,  it  does  appear  that  the  rotor 
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would  eventually  become  stable.  As  the  external  damping  increases,  the  speed  at  which 
the  system  becomes  stable  approaches  the  origin.  One  conclusion  from  Figure  23  is  that 
even  though  the  rotor  internal  damping  is  speed  independent  (due  to  the  internal 
hysteretic  damping),  the  stability  of  the  total  system  is  indeed  speed  dependent  due  to  the 
external  viscous  damped  bearings.  Figure  22  is  consistent  with  the  conclusions  reached 
by  Lund,  [11]  and  Dimentberg,  [15],  but  Figure  23  gives  some  new  results  on  the 
interaction  between  external  viscous  damping  and  internal  hysteretic  damping. 


Figure  23:  Log  Decrement  versus  Spin  Speed  for  System  with  Internal  Hysteretic 
Damping  and  a  Family  of  Constant  External  Viscous  Damping  Curves 

The  final  configuration  involved  adding  hysteretic  external  damping  to  a  rotor 
system  with  only  hysteretic  internal  damping.  Although,  the  FEM  model  was  not 
modified  to  represent  the  external  hysteretic  damping,  it  could  be  argued  that  one  could 
draw  the  plot  based  on  intuitive  reasoning.  First,  since  both  external  and  internal 
damping  are  hysteretic,  or  speed  independent,  the  logarithm  decrement  of  the  eigenvalue 
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solution  versus  the  spin  speed  would  be  horizontal  lines.  (Note:  In  the  next  section,  the 
FEM  does  actually  support  this  claim  for  the  case  of  zero  external  hysteretic  damping. 

See  Figure  31).  Therefore,  similar  to  the  previous  figure,  Figure  24  plots  the  log 
decrement  of  the  eigenvalue  solution  for  fixed  internal  hysteretic  damping  as  a  function 
rotor  spin  speed  Q.  But  this  time,  the  family  of  curves  represents  the  constant  external 
hysteretic  damping.  Note  that  as  the  external  hysteretic  damping  is  increased,  the  parallel 
lines  shift  toward  the  origin  and  toward  stability. 

Unlike  the  external  viscous  damping  where  the  spin  speed  determines  stability  (no 
matter  how  much  external  damping  is  applied),  the  external  hysteretic  damping  is  the 
only  factor  in  determining  stability  for  this  case.  Figure  24  is  consistent  with  the 
conclusions  reached  by  Lund,  [11]  and  Dimentberg,  [15],  although  it  clarifies  them  in  the 
respect  to  the  influence  of  external  hysteretic  damping  to  the  stability  of  the  system. 
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Figure  24:  Log  Decrement  versus  Spin  Speed  for  System  with  Internal  Hysteretic 
Damping  and  a  Family  of  Constant  External  Hysteretic  Damping  Curves 
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One  can  now  draw  the  following  conclusions: 


-  If  the  external  and  internal  damping  are  both  viscous  in  nature,  then  the  rotor  will 
always  be  stable  in  its  sub-critical  speed  range.  Stability  can  be  extended  into  the 
supercritical  regime  by  adding  external  damping  according  to  the  results  obtained  by 
Smith  in  equation  (56). 

-  If  the  external  and  internal  damping  are  both  hysteretic,  the  rotor  will  be  unstable  for 
all  speeds  unless  the  external  damping  is  large  enough  to  establish  stability.  In  this 
case,  the  rotor  will  be  stable  for  all  speeds. 

-  If  the  external  damping  is  viscous  and  the  internal  damping  is  hysteretic,  the  rotor 
will  be  unstable  in  a  range  of  speeds  beginning  at  £2  =  0+.  As  the  external  viscous 
damping  increases,  this  speed  range  shrinks  toward  £2  =  0  until  the  external  damping 
reaches  a  value  that  results  in  stability  for  all  speeds,  i.e.  for  £2  >  0. 
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Chapter  4  Finite  Element  Model  of  Rotor 


4.1  Finite  Element  Method  -  A  Timoshenko  Rotor  Element 

In  recent  years,  the  finite  element  method  has  been  successfully  applied  to 
rotordynamic  systems.  Although  all  rotor  elements  are  based  on  the  simple  beam 
element,  a  number  of  variations  have  been  developed.  One  such  element  is  based  on 
Timoshenko  beam  theory.  This  element  has  2  nodes  and  8  Degrees  of  Freedom  and  uses 
third  order  shape  functions  to  describe  the  bending  of  the  elements.  All  additional  inertia 
are  assumed  to  be  rigid  discs  with  lumped  mass  properties  and  the  bearings  are  assumed 
to  be  linear  and  discrete.  The  model  includes  rotary  inertia,  gyroscopic  moments,  and 
shear  deformation  effects.  A  more  detailed  discussion  of  a  Timoshenko  rotor  element  is 
located  in  Appendix  B. 

In  the  late  seventies,  Nelson  published  several  papers,  [16,  17],  to  determine  the 
accuracy  of  the  Timoshenko  rotor  element  and  to  document  current  works.  His  technique 
used  the  Lagrangian  approach  and  was  based  on  the  potential  and  kinetic  energy  of  the 
rotating  element  in  Figure  25. 
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X 


Figure  25:  Rotor  System  Configuration,  Nelson 

The  potential  energy  of  the  rotor  was  stored  in  two  forms,  bending  deformation,  and 
shear  deformation  (axial  loads  were  neglected)  and  was  given  by: 


where  the  shear  deformation  factor  is  given  by  K  and  the  second  derivative  of  the  bending 
deformation  in  the  y  direction  is  represented  by  V”end .  The  rotor  kinetic  energy  included 

Timoshenko  effects  of  rotary  inertia,  shear  deformation,  and  the  gyroscopic  energy.  It 
was  given  by: 


/  ,  i 

+  J-  pi  pa2ds  -  J  lpCiptdz 

0  2  o 


Evaluating  these  integrals  leads  to  the  familiar  forms  for  equations  (64)  and  (65): 
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(67) 


t=^J  ?}-nfe7  tc])fc}+|/,aJ 

Nelson  then  proceeded  to  systematically  derive  the  element  matrices  in  equations  (66) 
and  (67)  (see  Appendix  B).  By  the  application  of  Hamilton’s  extended  principle,  and 
using  the  above  equations,  Nelson  produced  the  following  undamped  matrix  equations  of 
motion: 

([MlraJHMrJ){q}-Q[G]{q}HCM}  +  [K]{q}^{f}  (68) 

where: 

{q}  -  fixed  frame  physical  coordinates 
{f}  -  fixed  frame  external  forces 
[Ce]  -  external  damping  matrix 

Similar  to  the  Jeffcott  rotor  with  the  non-central  disc,  the  eigenvalues  for  the 
undamped  equation  of  motion  (where  [Ce]  =  0)  occur  in  conjugate  pairs.  The  eigenvalues 
are  given  by:  Sj  =  ±iO)j ,  where  the  imaginary  pair  ±co  represents  the  forward  and 

backward  whirl  frequencies  of  the  shaft. 

4.2  Addition  of  Internal  Damping 

Based  on  the  work  of  Zorzi  and  Nelson,  [2],  the  previous  derivation  of  the 
rotordynamic  equations  of  motion  can  be  extended  to  include  the  contributions  of  internal 
damping.  Nelson  used  both  of  the  internal  viscous  and  hysteretic  damping  models 
discussed  in  the  previous  chapters.  Adding  the  constitutive  relationship  of  equations  (19) 
and  (36)  yields  the  following  equation: 
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cr=E 


(69) 


After  substituting  the  modified  constitutive  relationship  into  the  internal  bending 
moments  and  performing  the  integrations,  the  moments  can  be  expressed  as: 


where  the  internal  damping  matrix  was  defined  as: 


Now  placing  the  energy  contributions  from  these  moment  equations  into  the 
appropriate  kinetic  or  potential  energy  equations  (Lagrangian  approach),  Nelson  finally 
proceeds  to  derive  the  equation  of  motion  for  a  Timoshenko  shaft  finite  element  with 
internal  damping.  This  is  given  by: 

Wtrans ]  +  [M „,]){$}  +  (7/v  [K]  -  Q[G]  +  [Ce  ])  {q }  + 

^[K]+{-nva+-^=)[K,lr]{q)  =  {f)  (72) 

■Ji+%  v+n 

All  of  the  instabilities  of  the  system  are  characterized  in  the  new  skew-symmetric 
circulation  stiffness  matrix,  [2].  It  is  given  by: 
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skew 


sym 


0 


-12 

61 
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0 
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0 

0 

61 

-12 

0 

0 

61 


0 

-4/ 2  0 

61  0 
0  61 
0  2/ 2 

-2 12  0 


0 

-12  0 
-61  0  0 
0  -61  -4/ 2  0 


(73) 


(Incidentally,  the  circulation  matrix  given  in  reference  [2]  was  incorrect  due  to  misplaced 
negative  signs.  Equation  (73)  is  the  correct  circulation  matrix.)  It  is  also  important  to 
point  out  the  significance  of  therjv[K]{q]  term  in  equation  (72).  Although  this  term 

involves  the  viscous  internal  damping  coefficient,  it  does  not  produce  instabilities.  Recall 
Figure  16  where  the  direction  of  the  internal  viscous  damping  force  depended  on  the  sign 
of  the  (Q-(0)  term  in  equation  (53).  Equation  (72)  takes  this  effect  into  account  by  the 
7jv[K]{q}  term.  At  spin  speeds  less  than  the  critical  speed  (for  systems  without  external 

damping),  this  term  dominates,  resulting  in  a  stable  system.  But  at  speeds  above  the 
critical  speed,  the  terms  involved  with  the  circulation  matrix  dominate,  resulting  in  an 
unstable  system.  This  is  in  agreement  with  the  findings  documented  by  Smith  [6]  and 
Ehrich  [13],  and  is  also  consistent  with  Figure  16  in  Chapter  3. 

Solving  the  equations  of  motion  with  internal  damping  results  in  eigenvalues  in 
the  now  familiar  form  of:  Sj  =  A .  +  iCGj ,  where  again,  co  provides  shaft  whirl  frequency 


and  X  provides  the  orbit  growth/decay  rate  after  a  perturbation.  Nelson  then  followed 
Lund’s  logarithm  decrement  approach  in  determining  where  the  rotor  system  became 
unstable.  (Recall  that  the  instability  threshold  occurs  when  the  5  crosses  zero  when 
approaching  from  positive  values.)  The  Campbell  Diagram  in  Figure  26  shows  a  rotor 
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with  both  internal  viscous  and  hysteretic  damping  and  external  damping.  It  gives  a  good 
indication  of  how  the  stability  of  the  rotor  changes  as  the  values  of  8  approach  and 
surpass  zero  (at  the  instability  threshold). 


Rotor  Spin  Speed,  D. 


Figure  26:  Campbell  Diagram  for  Rotor  with  Internal  and  External  Damping 


4.3  Validation  of  FEM  Model 

In  line  with  its  work  on  magnetic  bearings,  The  Charles  Stark  Draper  Laboratory 
has  produced  a  rotor  dynamic  finite  element  model  (Antkowiak)  that  is  based  on 
undamped  Timoshenko  beam  theory.  In  this  model,  all  additional  inertia  are  assumed  to 
be  rigid  discs  with  lumped  mass  properties  and  the  bearings  are  assumed  to  be  discrete, 
undamped,  and  linear.  The  model  includes  rotary  inertia,  gyroscopic  moments,  and  shear 
deformation  effects. 
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After  making  the  appropriate  modifications  to  include  internal  and  external 
damping,  several  test  cases  were  run  to  check  the  accuracy  of  the  extended  rotor  model. 
The  bulk  of  the  test  cases  were  taken  directly  from  Nelson,  [2].  The  following  gives 
details  on  some  of  the  validation  of  the  extended  rotor  FEM  model. 

The  first  case  involved  an  undamped  (no  internal  or  external  damping)  steel  rotor 
supported  by  identical  bearings.  The  physical  dimensions  of  the  rotor  were  a  length  of 
1.27  meters  and  a  diameter  of  10.16  cm.  The  stiffness  of  the  undamped  isotropic 
bearings  was  1.75xl07  N/m. 

Figures  27  and  28  show  the  agreement  between  the  Draper  extended  model  and 
Nelson’s  solution.  The  first  critical  speed  given  by  Nelson  was  4950  rpm,  while  the 
Draper  FEM  model  calculated  4976  rpm.  The  small  discrepancies  could  be  attributed  to 
the  fact  that  this  particular  example  by  Nelson  did  not  take  into  account  shear  effects. 


Figure  27:  Draper  FEM  Model  Compared  to  Nelson  Solution, 
Undamped  Rotor  Supported  by  Identical  Bearings,  Modes  1  and  2 
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Figure  28:  Draper  FEM  Model  Compared  to  Nelson  Solution, 
Undamped  Rotor  Supported  by  Identical  Bearings,  Modes  1  and  2  (Detailed  View) 


The  next  test  case  involved  the  same  rotor  but  also  included  internal  viscous 
damping  in  the  rotor  (no  external  damping).  The  amount  of  internal  viscous  damping 
was  r|v  =  0.0002  sec'1.  Figure  29  shows  the  logarithm  decrements  of  the  first  and  second 
modes  (both  forward  and  backward)  as  a  function  of  spin  speed.  The  positive  values  of 
the  log  decrement  indicate  that  the  system  is  stable,  while  the  negative  values  indicate 
instability.  Good  agreement  is  shown  between  the  Draper  FEM  extended  model  and 
Nelson. 


Figure  29:  Log  Decrements,  Rotor  with  Viscous  Internal  Damping, 
No  External  Damping,  Modes  1  and  2 


4-8 


It  is  important  to  point  out  that  the  log  decrements  for  the  all  of  the  backward 
modes  are  increasing  and  therefore  will  always  be  stable.  On  the  other  hand,  the  forward 
modes  initially  are  stable  for  speeds  under  their  respective  critical  speeds,  then  become 
unstable  as  the  log  decrement  crosses  the  x-axis.  Note  that  this  agrees  with  the  theory  of 
internal  viscous  damping  described  in  Chapter  3. 

Figure  30  shows  the  logarithm  decrements  of  the  first  and  second  modes  when 
external  viscous  damping  is  added  to  the  rotor  system  through  the  bearings.  The  exact 
amount  of  external  damping  added  to  the  bearings  was  1.75x10  Ns/m.  Close  agreement 
is  again  shown.  Note  that  the  stability  of  the  system  was  improved,  i.e.  the  instability 
threshold  was  moved  approximately  5,000  rpm  higher  than  the  first  critical  speed.  In 
fact,  the  second  forward  mode  will  remain  stable  past  15,000  rpm.  This  effect  of  adding 
external  damping  to  the  system  also  agrees  with  the  internal  viscous  damping  model  in 
Chapter  3  and  in  Figure  21. 


Figure  30:  Log  Decrements,  Rotor  with  Viscous  Internal  Damping  and 
Viscous  External  Damping 


The  next  validation  involved  returning  to  the  same  rotor  and  this  time  only 
included  internal  hysteretic  damping  (T)h  =  0.0002).  Figure  31  shows  good  agreement 
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between  the  log  decrements  calculated  by  the  Draper  FEM  extended  model  and  Nelson. 
Similar  to  the  results  calculated  by  Lund  and  Dimentberg  (described  in  Chapter  3),  all  of 
the  forward  modes  are  unstable  and  all  of  the  backward  modes  are  stable. 


X  10* 

X  10* 

0.8 

Mode  1 

8 

Mode  2 

6 

0.6 

. o _ e — o — o — ©— a — o — o — © — o— -o — © — o — 0 — o  ' 

4 

0.4 

■ 

1  0-2 
E 

i  o 

STABLE 

c  2 
5 

E 

E  n 

STABLE 

3 

UNSTABLE 

1° 

£-0.2 

* 

UNSTABLE 

•0.4 

• 

*4 

x  -  Nelson  Forward  Mode 

o  -  Nelson  Backward  Mode 

-0.6 

- 

x  -  Nelson  Forward  Mode 

-0.8 

-1 

o  -  Nelson  Backward  Mode 

-6 

■ii - * - ■ - '  »i _ , _ , _ i 

0  5000  10000  15000  0  5000  10000  15000 

Rotor  Spin  Speed  (rpm)  Ro}or  Spln  Sp0Bd  (rpm) 


Figure  31:  Log  Decrements,  Rotor  with  Hysteretic  Internal  Damping, 
No  External  Damping,  Modes  1  and  2 


Figure  32  shows  the  logarithm  decrements  of  the  first  and  second  modes  when 
external  viscous  damping  is  added  through  the  bearings  to  this  rotor  system.  This 
addition  of  1.75xl03  Ns/m  per  bearing  completely  overcame  the  effects  of  the  internal 
hysteretic  damping  (r|H  =  0.0002),  resulting  in  a  rotor  system  stable  for  all  speeds  (both 
forward  and  backward  modes). 
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Figure  32:  Log  Decrements,  Rotor  with  Hysteretic  Internal  Damping 
and  Viscous  External  Damping 
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The  extended  rotor  FEM  model  was  successfully  validated  with  the  test  cases 
provided  by  Nelson,  [2].  Now  it  could  be  used  to  predict  the  critical  speeds  and 
instability  thresholds  for  other  rotors. 
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Chapter  5  Magnetic  Bearing  Simulation  Design 


5.1  State-Space  Model  of  Rotor  -  No  Internal  Damping 

In  addition  to  calculating  the  critical  speeds,  the  original  Draper  Lab  FEM 
program  (Antkowiak)  was  also  designed  to  produce  state  space  matrices  that  described 
the  rotor  system.  This  state  space  model  was  then  placed  in  a  Matlab™  magnetic  bearing 
simulation  developed  by  Draper  Laboratory  (Scholten)  to  test  the  closed  loop  control  of 
the  entire  rotor  system.  Originally,  the  program  produced  a  state-space  model  that  did  not 
include  internal  damping  and  only  included  a  small  amount  of  external  damping 
estimated  by  the  relationship: 

[Ce  ]  ~  2[£®] ,  diagonal  (74) 

where  the  external  viscous  damping  ratio  was  given  by 

The  original  FEM  model  used  generalized  eigenvectors,  modal  coordinates  and 
state  space  coordinates  to  calculate  the  state-space  model.  First,  the  undamped,  non¬ 
rotating,  homogeneous  equation  of  motion  in  physical  coordinates, 

[M]{$} +[*]{*}  =  [<>]  (75) 

was  used  to  find  a  set  of  eigenvectors  and  eigenvalues.  These  eigenvectors  were  then 
placed  in  then  xn  modal  matrix  [<&],  where  each  column  represented  a  different 
eigenvector  or  mode. 
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(76) 


pO)  p(2)  ...  p(«) 


The  eigenvectors  were  then  mass  normalized  by  the  following  equation: 


{Q(1)}  {Q(2)} 

V(Oa>  }r  [Af  Io0> } 


{Q(/i)} 


(77) 


If  one  converts  the  entire  non-homogeneous  differential  equation  derived  by  Nelson, 


m,raJ  +  [Mrol]){q}-Q[G]{q}+[Ce]{q}  +  [K]{q}  =  {f}  (69) 


to  modal  coordinates,  it  becomes: 

([Mtram  ]  +  [M  rot  ])[0]  {£}  +  (iCJ  -  OLGim  {£}  +  [K]  [O]  {£}  =  {/(0 }  (78) 

where  the  output  in  physical  coordinates  q  is  given  by  modal  coordinates,  in  the  form 
of: 


{q}  =  m{£}  (79) 

Now,  if  equation  (78)  is  pre-multiplied  by  the  transpose  of  the  mass  normalized  modal 
matrix,  [Of ,  it  gives 

[®]r([M,r.,J+tM„])[®]{|')+[®]r(iC,]-f!lG])[*]{^)  + 

+[*fm[2](f)=[®fi/(0) 

This  form  of  the  equation  lends  itself  to  the  use  the  property  of  orthogonality  where: 

[$f  + W,„  ll®]  =  [/] ,  diagonal 

[$f  (-£2[G]  +  [C,  ])[$]  =  ([G + 2[j>]])  (81  a,b,c) 

[Of  [^][0]  =  [t»2],  diagonal 
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At  this  point,  a  modal  reduction  technique  was  used  to  reduce  the  number  of  states 
in  the  mass,  stiffness,  and  gyroscopic  matrices  of  equations  (81).  This  eliminated  the 
higher  frequency  modes  outside  the  speed  range  of  interest. 

When  the  property  of  orthogonality  is  used,  equation  (80)  reduces  to: 

{f }  +  (2  ;[£»]  +  [G]){£ }  +  [O)2  ]  {£}  =  mT  [L]{f(t)}  (82) 

where  again,  the  external  damping  is  estimated  by: 

[CJ»2[£a],  diagonal  (74) 

It  is  important  to  emphasize  that  this  is  an  approximation  since  [<&]r[CJ[0]  will  not,  in 
general,  be  diagonal. 

If  one  lets  m  forces  be  applied  to  rotor  at  the  appropriate  nodal  points,  the  force 

vector  force  vector  { f(t) }  in  equation  (82)  must  be  modified  to  account  for  the  difference 

between  the  number  of  forces  applied,  m ,  and  the  number  of  modes  in  the  modal  matrix, 

n .  This  inflation  is  accomplished  by  then  x  rh matrix  [L]  where 

[L]  {/}  - >  [L){f}  (g3) 

nxm  m x  1  nxl 

To  construct  the  state  space  model,  the  equations  of  motion  must  now  be 
converted  to  state  space  coordinates  of  the  form: 

{ct}=m,{£}f  (84) 

With  this  substitution,  the  equation  of  motion  (82),  in  matrix  notation,  becomes: 


[[0]  [/]]{«:}  +  [[t»2]  [(2[Co)]  +  [G])\a}=mT[L]{f(t)} 

If  equation  (85)  is  combined  with  the  identity 

[[/]  [0]]{d}  +  [[0]  ~[I)]{a}=  0  (86) 
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it  gives: 


[/] 

[0] 


[0] 

[/] 


{a}  + 


[0]  -[/] 

[ co2]  -(2 [£a>]+\G]) 


{cc}=mT[L]{f(t)} 


(87) 


Now  using  the  standard  definitions  for  state  space  theory,  equation  (87)  can  also  be 
written  as 


{a)  =  [A]{a}+[B]{f(t)} 


(88) 


where  the  state  space  dynamics  matrix  [A]  is  defined  as 


[Ay- 


10]  [I] 

-[a2]  -(2[£y]+[G]) 


,2nx2n 


and  the  state  space  actuator  matrix,  [B],  is  defined  as 


[B]  = 


[0] 

mT[L] 


,2  nxm 


(89) 


(90) 


To  recover  the  physical  displacements  of  the  output  at  all  nodes  (or  locations  defined  by 
the  actuator  matrix),  state  space  theory  requires  that 

[q)  =  [C]{a]  (91) 

where  the  sensor  matrix,  [C],  is  defined  as: 

[C]  =  fc®]  [0]]  mx2n  (92) 


Although  this  definition  of  the  sensor  matrix  requires  the  displacement  sensors  to  be  co¬ 
located  at  the  positions  that  the  forces  are  applied,  the  original  model  and  simulation  did 
take  into  account  the  fact  that  the  sensors  were  not  exactly  collocated  with  the  magnetic 
bearing  actuators. 


5-4 


5.2  State-Space  Model  of  Rotor  -  Internal  Damping  Included 


If  internal  damping  is  included  in  the  FEM  model  and  equations  of  motion 
(creating  the  new  extended  rotor  model),  the  method  to  calculate  the  state-space  model 
remains  the  same,  but  the  results  are  more  complicated.  Starting  with  the  same 
undamped,  non-rotating,  homogeneous  equation  of  motion  in  physical  coordinates, 

[M]{q}  +  [K]{q}  =  [0]  ,  (75) 


a  set  of  eigenvalues  and  eigenvectors  are  calculated.  These  eigenvectors  were  then  placed 
in  the  modal  matrix  [O],  and  were  mass  normalized  by  equation  (77).  If  one  converts  the 
entire  non-homogeneous  differential  equation  derived  by  Nelson, 


m„.„  ]+[«„,]){</)+ <7)„m-Q[G] + [c,  ]>{«)+ 

{?>={/«)) 


1+''"  lK]+^vCl+-i3s=)lKdr] 


+  V  H 


fi+v F 


to  modal  coordinates,  it  becomes: 


([Mtmns  ]  +  [Mrot  ])[£]  {£} }  +  (: rjv  [K]  -  Q[G]  +  [Ce  ])[0]  {|}  + 

[&){£}  =  {/(»} 


l  +  V“  m  +  (77vQ+-^)[*c,] 


■\!^+vh 


(71) 


(93) 


where  the  output  in  physical  coordinates,  q,  is  given  by  modal  coordinates,  in  the  form 


of: 


{<?}  =  [£]{£} 


(79) 


Now,  if  equation  (93)  is  pre-multiplied  by  the  transpose  of  the  mass  normalized  modal 
matrix,  [<|>]r ,  it  gives 
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(94) 


mT  Wtrans  1  +  WM  ])[*]  {£}  +  [&?  fry  [AT]  -  Q[G]  +  [C,  ])[®]  {£}  + 


1+77// 
V1  +  /7« 


[AT]  +  (t7vQ  + 


77// 


V? 


+  77h 


®{|}  =  [^f{/(7)} 


This  form  of  the  equation  lends  itself  to  the  use  the  property  of  orthogonality  where: 


mT  mtrans  ]  +  [Mrot  ])[0]  =  [/],  diagonal 

mT  (77v  [K]  ~  Q[G)  +  [Ce])m  =  [G] 


(95a, b,c) 


m 


1  +  Vh  :[K]  +  (7jvQ+-fiM=)[Kcir] 


Vl  +  77// 


>/1  +  77// 


[<>]  =  [£] 


At  this  point,  a  modal  reduction  technique  was  used  to  reduce  the  number  of  states 
in  the  mass,  stiffness,  and  gyroscopic  matrices  of  equations  (95).  This  eliminated  the 
higher  frequency  modes  outside  the  speed  range  of  interest. 

When  the  property  of  orthogonality  is  used,  equation  (94)  reduces  to: 

{%}  +  ([GMt}  +  [Km  =  mT[L]{f(t)},  (96) 


where  this  time,  the  external  damping  is  not  approximated  by  the  damping  ratio.  Also, 
note  that  the  new  [K]  matrix  does  not  have  the  advantage  of  being  diagonal  (which 
complicates  the  solution). 

Again,  if  one  lets  m  forces  be  applied  to  the  rotor  at  the  appropriate  nodal  points, 
the  force  vector  force  vector  { f(t) }  in  equation  (96)  must  be  modified  to  account  for  the 
difference  between  the  number  of  forces  applied,  m ,  and  the  number  of  modes  in  the 
modal  matrix,  h .  This  inflation  is  accomplished  by  then  x  m matrix  [L]  where 


m  {/)  — *  [W}  (83) 

nxm  mx  1  nxl 

To  construct  the  state  space  model,  the  equations  of  motion  must  be  converted  to 
state  space  coordinates  of  the  form: 


5-6 


<«)=[{!),  (Ilf 


(84) 


To  recover  the  physical  displacements  of  the  output  at  the  nodes,  state  space 
theory  requires  that 

{■ q}  =  [C]{a }  (91) 

If  one  assumes  that  displacement  sensors  where  located  exactly  where  the  forces  where 


applied,  the  sensor  matrix  [C]  is  defined  as: 
[C]  =  [[0]  [0]]mx2n 
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(92) 


Note  that  the  [B]  and  [C]  matrices  remain  the  same  if  internal  damping  is  added  to  the 
system.  Only  the  state  space  dynamics  matrix,  [A],  is  affected. 

One  final  item  to  discuss  is  the  speed  dependence  of  the  [A]  matrix  for  either  case 
(with  or  without  internal  damping).  This  means  that  a  unique  [A]  matrix  exists  for  every 
spin  speed  Q.  To  overcome  this,  a  fitting  function  based  on  the  quadratic  Lagrange 
Polynomial  was  used  in  the  improved  simulation  (originally  by  Scholten)  to  calculate  the 
individual  [A]  matrices  as  the  speed  varied. 

5.3  Magnetic  Bearing  Simulation  Description 

After  the  FEM  program  calculated  the  state  space  matrices,  they  were  placed  in  a 
Matlab™  magnetic  bearing  simulation  developed  by  Draper  Laboratory  (Scholten)  to  test 
the  closed  loop  control  of  the  entire  rotor  system.  It  used  numerical  integration  to 
produce  time  history  plots  of  rotor  displacement  in  mils,  bearing  forces  in  lb,  bearing 
slew  rates  in  lb/s,  etc.  The  numerical  integration  was  performed  using  a  Runge-Kutta 
method  with  a  variable  step  size  ranging  from  1  x  10~8  to  1  x  10‘3  seconds. 

The  original  simulation  included  the  state  space  model  for  a  rotor  without  internal 
damping,  mathematical  models  for  the  bearings,  disturbance  functions,  and  rotor 
imbalance  and  included  appropriate  time  delays/lags  to  capture  the  workings  of  the  entire 
system.  Figure  33  shows  a  simplified  block  diagram  based  on  the  Simulink™  model 
(Scholten). 

The  flux-feedback  shown  in  the  figure  attempts  to  compensate  for  the  magnetic 
materials  non-linearities.  The  flux-feedback,  actuator  magnetics  and  power  amplifiers  are 
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assumed  to  work  together  as  a  linear  actuator,  modeled  as  having  force  and  slew  rate 
limits  combined  with  a  low  pass  frequency  response,  [18]. 

Due  to  the  modal  reduction  technique  mentioned  in  the  previous  section,  the  state 
space  matrices  used  by  the  simulations  only  included  frequencies  up  to  the  first  critical 
speed.  By  eliminating  the  higher  frequency  modes  outside  the  speed  range  of  interest, 
simulation  calculations  were  completed  faster. 


Actual 

Rotor 

Position 


Figure  33:  Simplified  Block  Diagram  of  Magnetic  Bearing  Simulation,  (Scholten) 
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5.4  Controller  Design 


The  bearing  controllers  used  for  the  simulation  were  designed  in  the  frequency 
domain  using  linear  Matlab™  based  software  tools  developed  at  Draper  Laboratory.  The 
control  system  individually  controlled  each  axis  on  the  bearing,  or  employed  a 
decentralized  control  strategy.  The  algorithm  used  for  the  controller  was  classical  PID. 
This  included  the  proportional  control  for  broad  band  stiffness,  the  integrator  control 
(with  anti-windup)  for  high  load  carrying  capacity,  and  derivative  control  to  damp 
disturbances,  [18]. 

The  actual  values  used  for  the  proportional  and  derivative  constants  were 
extremely  important  in  characterizing  the  maximum  stiffness  and  damping  characteristics 
of  the  bearings.  The  proportional  constant,  in  units  of  lb/mil,  provided  the  maximum 
stiffness  that  the  magnetic  bearings  could  provide.  The  derivative  constant,  in  units  of  lb- 
sec/mil,  gave  the  maximum  damping  that  the  bearings  could  provide. 
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Chapter  6  Modeling  and  Simulation  Predictions  and  Comparison 
to  Actual  Rotor  Data  -  No  Internal  Damping 

In  line  with  its  work  on  magnetic  bearings.  Draper  Laboratory  produced  a  rotor 
dynamic  finite  element  model  (Antkowiak),  an  accompanying  magnetic  bearing 
simulation  (Scholten),  and  a  rotor/magnetic  bearing  testing  apparatus.  To  determine  the 
accuracy  of  the  FEM  modeling  and  the  simulation,  several  test  runs  were  made  on  the 
rotor  test  apparatus.  Appendix  A  details  the  data  acquisition  rate  and  other  important 
information  regarding  actual  rotor  testing. 


6.1  FEM  Modeling 


The  actual  rotor  modeled  by  the  original  FEM  program  included  a  3/8”  diameter, 
14.7”  long  stainless  steel  shaft,  two  magnetic  actuator  rotors  and  an  elastic  coupling 
attached  to  one  end.  Figure  34  shows  the  basic  configuration. 


Actuator 


Elastic  Steel  Shaft 

Coupling 


Figure  34:  FEM  Model  of  Entire  Rotor 


The  pertinent  physical  properties  of  the  shaft  material  included  a  Modulus  of 
Elasticity  of  29.5  x  106  psi,  a  Poisson’s  ratio  of  0.33,  and  a  density  of  7.27  x  10'4  lbf/in3. 
The  shaft  was  modeled  as  20  Timoshenko  rotor  elements. 

As  Figure  34  shows,  an  elastic  coupling  was  necessary  to  connected  the  shaft  to 
the  drive  motor.  This  elastic  coupling  was  modeled  as  the  first  element  in  the  FEM 
model  with  a  Modulus  of  Elasticity  of  10.0  x  106  psi,  a  Poisson’s  ratio  of  0.45,  and  a 
density  of  7.50  x  10'4  lbf/in3,  [19]. 

The  two  magnetic  actuator  rotors  consisted  of  a  magnetic  alloy  laminated  onto  a 
common  aluminum  housing  secured  by  a  brass  nut.  Each  was  modeled  geometrically  as 
a  stiffened  shaft  with  added  mass  and  gyroscopic  inertia.  This  meant  that  each  of  these 
rotor  sections  would  have  the  same  diameter  as  the  shaft  (in  the  FEM  input  file),  but 
would  have  an  artificially  increased  stiffness  to  compensate  for  the  effects  of  the  larger 
diameter  rotor.  This  increased  bending  stiffness  was  calculated  to  be  4.69  times  that  of 
steel.  Due  to  the  interference  fit,  only  half  of  this  value  was  used  to  modify  the  Modulus 
of  Elasticity.  This  resulted  in  a  Modulus  of  Elasticity  of  69.3  x  106  psi,  a  Poisson’s  ratio 
of  0.33,  and  a  density  of  8.40  x  10"4  lbf/in3,  for  the  elastic  coupling,  [19]. 

With  regard  to  the  additional  mass  and  inertia,  Figure  35  shows  the  locations 
where  they  were  concentrated  along  the  rotor  axial  length.  In  addition,  Table  2  shows  a 
summary  of  the  added  mass  and  inertia  for  the  magnetic  actuator  rotor.  The  input  file  for 
the  FEM  program  describing  this  rotor  is  located  in  Appendix  C. 
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Figure  35:  FEM  Model  of  Rotor  Actuator 


Table  2:  Calculations  of  Ad 

Ided  Mass  and  Inertia  Due  to  Actuator  Rotor 

Location 

1 

2 

3 

4 

Axial  Length  (1.3  in) 

0.18 

0.47 

0.22 

Mass  (188  lbm)xl0'b 

26. 

62. 

68. 

32. 

■■EjQnjljym 

1.2 

3.6 

4.1 

1.5 

Ixxx  lO^in4 

2.2 

5.3 

5.8 

2.7 

An  important  part  of  the  FEM  analysis  was  determining  the  stiffness  and  damping 
characteristics  of  the  bearings.  Although  magnetic  bearings  are  able  to  change  their 
stiffness  and  damping  depending  on  the  position  of  the  shaft,  the  FEM  program  required 
constant  values  for  its  calculations.  In  an  attempt  to  model  the  magnetic  bearings,  the 
stiffness  and  damping  was  estimated  to  be  50  lbs/in  and  1.5x1  O'2  lbs-sec/in,  respectively. 
These  were  based  upon  the  maximum  stiffness  and  damping  that  the  controller  could 
provide.  More  specifically,  they  were  taken  directly  from  the  proportional  and  derivative 
constants  (in  units  of  Ib/mil  and  lb-s/mil,  respectively)  of  the  PID  controller  with  the  time 
delays  and  system  lags  taken  into  account. 

The  rotor  free-free  natural  frequencies  are  presented  in  Table  3.  The  first  four 
modes  show  the  rigid  body  modes  at  zero  Hertz  since  both  rotor  ends  are  free.  (The 
additional  zero  frequency  modes  are  present  due  to  the  forward  and  backward 
components  of  the  state-space  model.)  These  rigid  body  modes  consist  of  the  conical 
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mode  and  the  cylindrical  (or  bounce)  mode.  The  last  entries  in  the  table  are  the  first  three 
natural  frequencies  of  the  non-rotating  rotor. 


Table  3:  Rotor  Free-Free  Natural  Frequencies  for  Q  =  0 


Mode# 

Rotor  Natural 
Frequencies  (Hz) 

1-2 

0.0 

3-4 

0.0 

5-6 

269.9 

7-8 

790.5 

9-10 

1,559.7 

Figure  36  presents  the  Campbell  Diagram  for  the  rotor  model.  The  symmetric 
bearings  have  a  stiffness  and  damping  of  50  lbs/in  and  1.5x1  O'2  lbs-sec/in,  respectively. 
The  intersections  between  the  diagonal  line  (synchronous  whirl)  and  the  whirl  frequency 
pairs  indicate  possible  critical  speeds.  The  rigid  body  modes,  the  whirl  frequency  pairs, 
and  the  critical  speeds  for  the  rotor  model  are  given  in  Table  4.  The  critical  speed  for 
each  mode  is  in  bold  text. 


0  0.5  1  1.5  2  2.5  3  3.5 

Rotor  Spin  Speed  (rpm)  x  104 

Figure  36:  Campbell  Diagram  for  Rotor 
Bearing  Stiffness  =  50  lbs/in,  Bearing  Damping  =  1.5xl0'2  lbs-sec/in 
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Table  4:  Rotor  Whirl  Frequency  Pairs 
Bearing  Stiffness  =  50  Ibs/in,  Bearing  Damping  =  1.5xl0'2  lbs-sec/in 


Mode 

Rotor  Whirl  Frequency  Pairs  (Hz) 

Rigid  Body 

42.1 

,  67.2* 

Backward 

Forward2 

1 

268.9 

272.5 

2 

779.7 

802.0 

3 

1,531 

1,590 

1 1 - 1 - 1 - 1 ^ - 

Rigid  body  modes  are  uncoupled,  Rotor  Critical  Speeds 


Note:  The  original  FEM  program  predicted  two  rigid  body  modes  (conical  and  bounce) 
and  gave  their  uncoupled  values.  In  reality,  the  two  rigid  body  modes  are  actually 
coupled.  This  coupling  typically  results  in  only  one  noticeable  rigid  body  mode  (a 
combination  of  both  conical  and  bounce  modes),  [19]. 

After  the  whirl  speed  analysis  was  complete,  the  original  FEM  program 
(Antkowiak)  generated  state  space  matrices  for  the  rotor  model.  The  bode  plot  of  the 
free-free  rotor  model  is  given  in  Figure  37.  It  shows  the  first  critical  speed  at 
approximately  1700  rad/sec  or  270  Hz,  which  is  consistent  with  Table  3. 


200 


io'1  io°  io'  io2  io3  io4 

Frequency  (rad/sec) 


Figure  37 :  Bode  Plot  of  Free-Free  Rotor  Model 
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Prior  to  entering  the  rotor  model  into  the  simulation,  the  root  loci  of  the  free-free 
rotor  model  was  plotted.  A  root  locus  plot  illustrates  how  the  closed  loop  poles  and  zeros 
of  a  system  change  as  the  gain  is  varied.  Although  proportional  gain  is  a  simple  choice 
for  a  controller,  it  gives  a  quick  look  at  the  system  stability,  gives  insight  on  possible 
controller  designs,  and  illustrates  the  complexity  of  the  dynamic  model. 

Figure  38  shows  the  how  the  root  loci  change  as  the  rotation  speed  changes.  Note 
that  due  to  the  scale  of  the  plots,  the  root  loci  lines  between  the  poles  (zero  gain)  and  the 
zeros  (infinite  gain)  are  not  visible.  At  zero  rotational  speed,  all  of  the  poles  and  zeros 
are  located  on  the  imaginary  axis.  As  the  rotation  speed  increases,  the  poles  and  zeros 
move  in  a  curved  path  toward  the  negative  real  axis.  Once  aligned  with  the  real  axis, 
some  begin  to  migrate  towards  negative  infinity  and  some  approach  (but  never  reach)  the 
origin.  Since  the  poles  of  this  system  never  cross  into  the  right  half  plane,  the  system 
should  be  stable  for  all  rotation  speeds.  In  addition,  since  the  zeros  of  the  system  never 
enter  the  right  half-plane,  it  indicates  that  a  simple  controller  might  be  sufficient  to 
control  the  system,  [20]. 
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Figure  38:  Root  Loci  of  Free-Free  Rotor  Model  as  Rotation  Speed  Increases 
(No  Internal  or  External  Damping) 


6.2  Magnetic  Bearing  Simulation  Results 


Using  the  calculated  state-space  matrices  from  the  original  FEM  model 
(Antkowiak),  the  magnetic  bearing  simulation  (Scholten)  predicted  the  critical  speeds  and 
rigid  body  modes  for  the  rotor.  Figure  39  shows  the  simulated  response  of  the  rotor  to  a 
0  to  325  Hz  ramp  in  one  second.  The  top  time  history  plot  shows  the  displacement  of  the 
rotor  at  bearing  one  (located  next  to  the  elastic  coupling)  in  both  the  vertical  and 
horizontal  axes.  The  second  time  history  shows  the  displacement  of  the  rotor  at  bearing 
two.  Underneath  the  time  histories,  additional  plots  are  provided  to  show  the  radial  shaft 
orbit  during  the  simulated  test  run.  The  dashed  lines  in  the  orbit  plots  represent  the 
touchdown  bearings  (8  mil  air  gap).  In  general,  the  results  of  Figure  39  indicate  that  the 
PID  controller  can  successfully  control  the  rotor  as  it  passes  through  the  first  critical 
speed  (occurring  at  approximately  0.92  seconds  or  290  Hz)  without  contacting  the 
touchdown  bearings.  This  also  confirms  the  assumption  based  on  the  root  loci  plots 
(Figure  38)  that  a  simple  controller  would  be  sufficient  to  control  this  system. 
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Figure  39:  Simulation  Prediction  of  Rotor  Response  during  a  0  -  325  Hz  Ramp 


To  get  a  better  insight  on  the  closed  loop  control  of  the  plant,  a  number  of  plots 
were  generated  using  frequency  domain  based  tools  developed  at  Draper  Laboratory. 
Figure  40  shows  the  frequency  response  of  the  rotor  located  at  bearing  number  two 
(farthest  away  from  the  elastic  coupling)  in  terms  of  magnitude  (output/input  ratio  in  dB) 
and  phase  (output  -  input  in  degrees)  for  a  rotor  speed  of  150  Hz.  This  spin  speed  was 
chosen  as  a  representative  operating  speed  for  the  rotor  due  to  motor  speed  limitations. 
The  first  response,  A,  shown  is  the  plant  model.  For  this  case,  the  plant  includes  bearing 
one  represented  by  a  first  order  lag  model.  The  second  response,  B,  shows  the  controller 
response,  which  includes  calculation  and  zero-order-hold  delays.  The  final  pair  of  plots, 
C,  shows  the  open  loop  (without  feedback)  response  of  the  plant  and  the  controller 
combined.  Phase  is  positive  from  approximately  40  to  300  Hz  indicating  the  frequency 
band  where  damping  occurs,  [22], 
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A:  PLANT,  C:\MATLAB4\BIN\P4S10.MAT 


Frequency  [Hz] 


Figure  40:  (A)  Plant  Response,  (B)  Controller  Response,  (C)  Plant  and 
Controller  Combined  (Open-Loop)  Response  for  a  Rotor  Spin  Speed  of  150  Hz 


Figure  41  is  a  Nichols  plot  of  the  magnitude  versus  the  phase  from  the  open-loop 
plant  and  controller  combined  case  shown  in  Figure  40.  The  rotor  speed  is  150  Hz.  The 
‘X’  indicates  the  operating  point  of  the  control  law  relative  to  the  complex  structure  of 
the  rotordynamics.  The  distance  of  the  ‘X’  from  the  line  gives  the  gain  margin  as  15  dB 
and  the  phase  margin  as  30  degrees.  These  adequate  margins  indicate  that  the  controller 
will  be  reliably  stable,  with  a  damped  response  to  disturbances,  [21]. 
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Figure  41:  Open  Loop  Gain  vs.  Phase  at  150  Hz 

Figure  42  presents  the  open  and  closed  loop  frequency  response  of  the  controlled 
bearing  furthest  away  from  the  motor  (bearing  2).  The  dashed  line  corresponds  to  the 
uncontrolled  rotor  (combined  plant  and  controller  response  of  Figures  40  and  41)  and  the 
solid  line  shows  the  controlled  rotor.  The  difference  between  the  two  represents  the 
action  of  the  magnetic  bearings  upon  the  rotor.  At  lower  frequencies  the  effect  is  large 
because  the  bearings  can  tightly  constrain  the  rotor.  At  high  frequencies  the  effect  is 
minimal  due  to  the  inherent  difficult  of  moving  the  rotor  since  the  displacement  is 
proportional  to  force/frequency2,  [21].  In  addition,  these  higher  frequencies  are  beyond 
the  actuator  bandwidth.  The  critical  speed  for  the  closed  loop  system  shown  in  this  figure 
is  289  Hz,  which  is  close  to  the  value  estimated  from  Figure  39. 
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Plant  and  Controller  Combined  Open/Clos ed-Loop  FRF:  [PID  kp=.05  ki=3.00  kd=0.002  pdhz=1000]  Gain  =  1 


Frequency  [Hz] 

Figure  42:  Open  and  Closed  Loop  Frequency  Response,  150  Hz 


To  find  the  rigid  body  modes  predicted  by  the  simulation,  a  Power  Spectral 
Density  plot  of  the  rotor  system,  at  a  constant  speed  of  148  Hz,  was  produced.  Similar  to 
a  Bode  plot,  a  PSD  takes  the  all  of  the  modes  acting  on  the  system  and  plots  them  as  a 
function  of  frequency.  Figure  43  shows  the  dominate  mode  of  148  Hz  (related  to  the 
rotor  spin  speed)  and  small  peak  occurring  at  289  Hz  (critical  speed).  Although  there  is 
not  a  distinct  peak,  the  relative  maximum  of  the  curve  gives  a  rigid  body  mode  estimate 
of  47  Hz.  A  comparison  of  the  FEM  model  predictions  to  the  closed  loop  simulation  is 
given  in  Table  5. 
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Figure  43:  PSD  of  Simulation  at  Constant  Rotor  Spin  Speed  of  148  Hz 


Table  5:  Comparison  of  Original  FEM  Model  to  Original  Simulation 
FEM  Model:  Bearing  Stiffness  =  50  Ibs/in,  Bearing  Damping  =  1.5xl0‘2  Ibs-sec/in 
Simulation:  Kp  =  0.05  lb/mil,  Ki  =  3.00  Ib/mil-s,  KD  =  0.002  lb-sec/mil, 


Derivative  Bandwidth  =  1C 

>00  Hz 

Mode 

FEM 

Prediction 

Simulation 

Rigid  Body  Mode 

42.1,67.2* 

47 

1st  Critical  Speed 

272.5 

289 

*  Rigid  body  modes  are  uncoupled. 


6.3  Actual  Rotor  Testing  and  Comparison  to  Simulation 


To  determine  the  accuracy  of  the  FEM  model  and  the  Matlab™  magnetic  bearing 
simulation,  several  test  runs  were  made  on  a  rotor  test  apparatus  located  at  Draper 
Laboratory.  Appendix  A  details  the  data  acquisition  rate  and  other  important  information 
regarding  actual  rotor  testing. 
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Figures  44  and  45  show  two  representative  simulation  predictions  compared  to 
actual  rotor  rig  test  data.  For  each  pair  of  plots,  the  simulation  predictions  are  located  on 
the  left  and  the  actual  rotor  test  runs  on  shown  on  the  right.  All  cases  were  run  at  a 
constant  rotation  speed  and  all  displacements  were  measured  in  mils  from  the  bearing 
centerline.  Figures  44  and  45  show  good  agreement  in  terms  of  the  predicted  and  actual 
displacements  for  bearing  number  2,  although  the  simulation  predicted  somewhat  lower 
displacements  for  the  bearing  near  the  elastic  coupling.  This  suggests  that  more  work 
needs  to  be  done  on  the  elastic  coupling  block  of  the  simulation.  In  terms  of  the  vertical 
displacements  at  the  bearings,  the  simulation  had  good  correlation  to  the  actual  rotor  data, 
while  the  predicted  horizontal  displacements  were  slightly  smaller  than  the  actual  rotor 
data. 


seconds 


Bearing  1- mils  Bearing  2 -mils 


Figure  44:  Predicted  and  Actual  Rotor  Displacements,  102  Hz 
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sim4.  Rotor  4.  148  Hz  Constant.  AB:afact=20.  0  lb  load.  24-Mar-99 


Actual  Short  Rotor  Data  - 148  Hz  03-Mar-l999 


Bearing  1  -  mils  Bearing  2  -  mils  Bearing  1  -  mils  Bearing  2  -  mils 


Figure  45:  Predicted  and  Actual  Rotor  Displacements,  148  Hz 


The  only  area  where  the  simulation  had  poor  accuracy  with  regard  to  the  actual 
rotor  displacements  was  in  the  speed  range  near  the  rigid  body  mode.  Since  the 
simulation  predicted  the  rigid  body  mode  at  47  Hz,  its  displacement  predictions  at  the 
actual  rigid  body  mode  of  52  Hz  would  not  be  as  large.  This  is  shown  in  Figure  46  with 
the  simulation  prediction  on  the  left  and  the  actual  rotor  data  on  the  right. 

slm4.  Rotor  4.  52  Hz  Constant.  AB:afact=20.  0  lb  load,  imbl  =0.001, 0.001  4-Mar-99  Actual  Short  Rotor  Data  -  52  Hz  03-Mar-1999 


Bearing  1 -mils  Bearing  2  -  mils  Bearing  1  -  mils  Bearing  2 -mils 

Figure  46:  Predicted  and  Actual  Rotor  Displacements,  52  Hz 
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A  more  important  measure  of  the  simulation  was  its  ability  to  predict  the  location 
of  the  rigid  body  modes  and  the  critical  speeds.  With  regard  to  the  first  critical  speed,  the 
actual  rotor  could  not  be  subjected  to  a  ramp  from  0  to  325  Hz  in  one  second  (like  the 
simulation).  In  addition,  hardware  limitations  prevented  the  rotor  from  achieving  speeds 
greater  than  200  Hz.  Therefore,  to  determine  the  first  critical  speed  and  get  another 
estimate  of  the  rigid  body  mode,  a  PSD  plot  was  made  from  actual  rotor  data  at  148  Hz. 
Figure  47  shows  the  dominate  mode  of  148  Hz  (related  to  the  rotor  spin  speed),  a  peak  at 
52  Hz  (rigid  body  mode)  and  a  small  peak  occurring  at  295  Hz  (critical  speed). 


Hertz 


Figure  47:  PSD  of  Actual  Rotor  at  Constant  Rotor  Spin  Speed  of  148  Hz 

Table  6  presents  a  summary  of  the  results  for  the  FEM  rotor  model,  the 
simulation,  and  the  actual  rotor  data.  The  FEM  critical  speed  prediction  was  slightly 
lower  than  the  actual  rotor  rig  critical  speed.  This  was  probably  due  to  the  simplicity  of 
the  FEM  model.  This  model  did  not  take  into  account  rotor  unbalance,  it  assumed  that 
the  bearing  and  properties  were  constant,  and  it  did  not  model  any  other  overall  system 
effects  such  as  magnetic  bearing  material  non-linearity,  actuator  lag,  etc.  (Note:  Since  the 
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FEM  produced  uncoupled  rigid  body  modes,  they  could  not  be  directly  compared  to  the 
actual  rotor  rig  data.) 


Table  6:  Comparison  of  Original  FEM,  Modified  Simulation  and  Actual  Rotor  Rig 
Data,  FEM  Model:  Bearing  Stiffness  =  50  Ibs/in,  Bearing  Damping  = 
1.5xl0'2  lbs-sec/in,  Simulation:  Kp  =  0.05  lb/mil,  Ki  =  3.00  lb/mil-s, 
_ Kp  =  0.002  lb-sec/mil,  Derivative  Bandwidth  =  1000  Hz _ 


Mode 

FEM 

Prediction 

Simulation 

Prediction 

Actual  Rotor 
Rig  Data 

42.1,67.2 

47 

52 

272.5 

289 

295* 

*  Indicates  condition  not  tested. 


In  contrast  to  the  simple  FEM,  the  high  fidelity  simulation  (originally  by  Scholten 
with  slight  modifications  described  in  Chapter  5)  predicted  a  first  critical  speed  that 
closely  matched  the  actual  rotor  test  data,  although  the  rigid  body  mode  predictions  were 
not  as  accurate.  With  regard  to  the  displacements  predicted  by  the  simulation,  good 
agreement  was  seen  in  terms  of  the  predicted  and  actual  displacements  for  bearing 
number  2,  although  the  modified  simulation  predicted  somewhat  lower  displacements  for 
the  bearing  near  the  elastic  coupling.  This  suggests  that  more  work  needs  to  be  done  on 
the  elastic  coupling  block  of  the  simulation.  In  terms  of  the  vertical  displacements  at  the 
bearings,  the  modified  simulation  had  good  correlation  to  the  actual  rotor  data,  while  the 
predicted  horizontal  displacements  were  slightly  smaller  than  the  actual  rotor  data.  The 
only  area  where  the  simulation  had  poor  accuracy  with  regard  to  the  actual  rotor 
displacements  was  in  the  speed  range  near  the  rigid  body  mode.  Since  the  simulation 
predicted  the  rigid  body  mode  at  47  Hz,  its  displacement  predictions  at  the  actual  rigid 
body  mode  of  52  Hz  would  not  be  as  large. 


6-17 


Chapter  7  Modeling  and  Simulation  Predictions  for  a  Rotor  with  Internal  Damping 

7.1  FEM  Modeling 

The  same  rotor  at  Draper  Laboratory  was  again  modeled  using  the  extended  FEM 
program.  But  for  these  cases,  a  number  of  different  combinations  of  internal  damping 
was  included  in  the  rotor.  First,  internal  viscous  damping  (r|v  =  2x1  O'6)  was  included. 
Next,  only  internal  hysteretic  damping  (tJh  =  2xl0'4)  was  included  in  the  rotor.  Finally, 
both  internal  hysteretic  and  viscous  damping  were  added  to  the  rotor. 

The  actual  amounts  of  internal  damping  that  were  included  in  the  rotor  were 
estimated  based  on  references,  [2,  10,  and  11].  Even  though  the  actual  amount  of  internal 
damping  in  the  rotor  remains  unknown,  these  estimated  values  are  still  useful  to  show 
how  the  simulation  (Scholten)  and  controller  would  react  to  the  destabilizing  effects  of 
internal  damping  in  general. 

An  important  part  of  the  FEM  analysis  was  determining  the  stiffness  and  damping 
characteristics  of  the  bearings.  Although  magnetic  bearings  were  able  to  change  their 
stiffness  and  damping  depending  on  the  position  of  the  shaft,  the  FEM  program  required 
constant  values  for  its  calculations.  In  an  attempt  to  model  the  bearings,  the  stiffness  and 
damping  was  assumed  to  be  50  lbs/in  and  1.5x10'  lbs-sec/in  respectively.  These  were 
based  upon  the  maximum  stiffness  and  damping  that  the  controller  could  provide.  More 
specifically,  they  were  taken  directly  from  the  proportional  and  derivative  constants  of 
the  PID  controller  (in  units  of  lb/mil  and  lb-s/mil,  respectively)  with  the  time  delays  and 
system  lags  taken  into  account.  The  FEM  calculated  whirl  frequency  pairs  and  estimates 
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of  the  instability  thresholds  for  the  various  combination  of  internal  damping  are  given  in 
Table  7.  The  critical  speed  for  the  first  mode  is  in  bold  text. 


Table  7:  FEM  Predicted  Rotor  Instability  Thresholds 
Bearing  Stiffness  =  50  lbs/in,  Bearing  Damping  =  1.5x1  O'2  Ibs-sec/in 


Internal  Damping 

FEM  Whirl  Frequencies 
(Hz) 

Instability  Threshold 
(Hz) 

Backward 

Forward 

T)v  =  2x1  O'6 

268.9 

441 

tjh  =  2xl04 

268.9 

Stable  for  all  Speeds 

rjv  =  2x1 0'6,  T|h  =  2x1 0‘4 

268.9 

272.5* 

425 

*  Critical  speed 


After  the  whirl  speed  analysis  was  complete,  the  FEM  program  generated  state 
space  matrices  for  the  three  rotor  models.  The  bode  plots  of  the  free-free  rotor  models 
(without  external  damping  or  stiffness)  are  given  in  Figures  48-50.  All  plots  show  the 
same  values  for  the  first  critical  speed  at  approximately  1700  rad/sec  or  270  Hz,  which  is 
consistent  with  Table  7.  In  addition,  all  plots  show  the  erratic  tendencies  of  the  phase  as 
the  frequency  is  increased.  This  could  be  a  result  of  adding  internal  damping  to  the 
system,  and  needs  further  investigation. 


Frequency  (rad/sec) 


Frequency  (rad/sec) 

Figure  48:  Bode  Plot  of  Free-Free  Rotor  Model  at  250  Hz 
with  Internal  Viscous  Damping  T|v  =  2xl0'6 
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Frequency  (rad/sec) 

Figure  49:  Bode  Plot  of  Free-Free  Rotor  Model  at  250  Hz 
with  Internal  Hysteretic  Damping  r|H  =  2xl0'4 


Frequency  (rad/sec) 


Figure  50:  Bode  Plot  of  Free-Free  Rotor  Model  at  250  Hz,  with 
Internal  Viscous  Damping  T|v  =  2xl0*6,  Internal  Hysteretic  Damping  riH  =  2x1  O'4 

Prior  to  entering  the  rotor  model,  into  the  modified  simulation,  the  root  loci  of  the 
free-free  rotor  models  were  plotted.  Figure  5 1  shows  the  how  the  root  loci  of  the  rotor 
with  internal  viscous  damping  change  as  the  rotation  speed  changes.  At  zero  rotational 
speed,  a  few  poles  and  zeros  are  located  at  the  origin,  with  the  rest  in  the  stable  left  half 
plane.  As  the  rotation  speed  increases,  the  two  pairs  of  poles  and  zeros  begin  to  migrate 
toward  the  imaginary  axes.  Between  250  and  300  Hz  the  pairs  actually  cross  over  into 
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the  unstable  region.  This  crossover  occurs  at  the  critical  speed  of  the  free-free  rotor. 
This  is  consistent  with  the  general  characteristics  of  internal  viscous  damping  (see  the 
discussion  in  Chapter  3).  Since  the  poles  and  zeros  of  this  rotor  do  cross  over  into  the 
right  half  plane,  this  indicates  that  a  rotor  with  internal  viscous  damping  might  be 
difficult  to  control  at  speeds  above  its  critical  speed,  [20]. 
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Figure  51:  Root  Loci  of  Free-Free  Rotor  Model  as  Rotation  Speed  Increases 
(Internal  Viscous  Damping  r|v  =  2xl0‘6) 

Figure  52  shows  the  how  the  root  loci  of  the  rotor  with  internal  hysteretic 
damping  change  as  the  rotation  speed  changes.  For  all  rotational  speeds,  a  pair  of  poles  is 
located  in  the  left  and  right  half  planes.  This  indicates  that  the  free-free  rotor  will  be 
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unstable  for  all  speeds.  These  results  are  consistent  with  the  general  characteristics  of 
internal  hysteretic  damping  (see  the  discussion  in  Chapter  3).  As  the  rotation  speed 
increases,  the  two  pairs  of  zeros  begin  to  migrate  toward  the  poles.  Since  this  system  has 
poles  and  zeros  located  in  the  right  half  plane,  this  indicates  that  a  rotor  with  internal 
hysteretic  damping  might  be  difficult  to  control  at  all  speeds,  [20]. 
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Rotor  Spin  Speed  =  500  Hz 


Figure  52:  Root  Loci  of  Free-Free  Rotor  Model  as  Rotation  Speed  Increases 
(Internal  Hysteretic  Damping  T)h  =  2xl0‘4) 


Figure  53  shows  the  how  the  root  loci  of  the  rotor  with  internal  viscous  and 
hysteretic  damping  change  as  the  rotation  speed  changes.  In  general,  the  trend  is  very 
similar  to  Figure  51  (internal  viscous  damping  only),  where  the  poles  and  zeros  migrate 
toward  and  cross  into  the  unstable  region.  The  difference  lies  in  that  the  system  can  be 
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unstable  at  all  speeds  depending  on  the  value  of  the  gain.  This  is  most  apparent  in  the  top 
left-hand  comer  of  Figure  53  (rotor  spin  speed  of  1  Hz).  Similar  to  the  internal  viscous 
damping  only  case,  these  results  indicate  that  a  rotor  with  both  internal  viscous  and 
hysteretic  damping  might  also  be  difficult  to  control  at  speeds  above  its  critical  speed, 
[20]. 

Although  it  may  appear  odd  that  this  free-free  rotor  (without  external  damping) 
can  be  stable  even  though  it  has  internal  hysteretic  damping,  one  must  recall  the 
importance  Figure  16  from  Chapter  3  or  the  stabilizing  term,  T}v[K]{q) ,  from  equation 
(72).  The  figure  shows  that  the  particular  direction  of  the  internal  viscous  damping  force 
depends  on  the  sign  of  the  (£2-co).  For  this  particular  rotor,  the  internal  hysteretic 
damping  was  in  the  same  direction  as  the  rotor  whirl  and  the  internal  viscous  damping 
was  in  the  opposite  direction.  So,  for  speeds  below  the  critical  speed  of  the  system,  the 
stabilizing  effects  of  internal  viscous  damping  counteracted  the  destabilizing  effects  of 
the  internal  hysteretic  damping. 
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7.2  Magnetic  Bearing  Simulation  Predictions  (Internal  Damping  Included) 


Using  the  calculated  state-space  matrices  from  the  extended  FEM  model,  the 
simulation  attempted  to  predict  the  instability  thresholds  for  the  various  combinations  of 
internal  damping.  Figure  54  shows  the  simulated  response  of  the  rotor  with  internal 
viscous  damping  during  a  0  to  400  Hz  ramp  in  one  second.  At  the  closed-loop  system’s 
critical  speed  of  289  Hz  (approximately  0.72  seconds)  the  rotor  becomes  unstable.  This 
instability  threshold  is  well  below  the  FEM  predicted  speed  of  441  Hz  found  in  Table  7. 


sim4i.  Rotor  4,  viscous  internal  damping,  0-400  Hz  Ramp.  25-Mar-99 
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Figure  54:  Simulation  Prediction  of  Rotor  Response  during  a  0  -  400  Hz  Ramp, 
Internal  Viscous  Damping  r|v  =  2x1  O'6 


To  get  a  better  insight  on  the  closed  loop  control  of  the  plant,  a  number  of  plots 
were  generated  using  frequency  domain  based  tools  developed  at  Draper  Laboratory. 
Figure  55  shows  the  frequency  response  of  the  rotor  located  at  bearing  number  two  in 
terms  of  magnitude  (output/input  ratio  in  dB)  and  phase  (output  -  input  in  degrees)  for  a 
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rotor  speed  of  250  Hz.  The  first  response.  A,  shown  is  the  plant  model.  For  this  case,  the 
plant  includes  bearing  one  represented  by  a  first  order  lag  model.  The  second  response, 
B,  shows  the  controller  response.  The  final  pair  of  plots,  C,  shows  the  open  loop  (without 
feedback)  response  of  the  plant  and  the  controller  combined.  Similar  to  the  bode  plot  in 
Figure  48,  note  the  sudden  phase  angle  changes  for  both  the  plant  and  for  the  combined 
open  loop  plant  and  controller. 
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Figure  55:  (A)  Plant  Response,  (B)  Controller  Response,  (C)  Plant  and 
Controller  Combined  (Open-Loop)  Response  for  a  Rotor  Spin  Speed  of  250  Hz, 
Internal  Viscous  Damping  T]v  =  2xl0'6 


A  number  of  different  PID  controllers  were  tested  in  an  attempt  to  stabilize  the 
rotor  above  its  critical  speed.  These  efforts  focused  on  controller  design  in  the  frequency 
domain  with  subsequent  testing  in  the  full  simulation.  Increasing  (or  decreasing)  the 
proportional,  integral,  and  derivative  constants  had  little  effect  on  the  stability  of  the 
system. 


7-10 


In  a  sense,  the  PDD  controller  was  “boxed  in”  by  the  complex  rotordynamics. 

This  is  shown  best  by  the  Nichols  plot  in  Figure  56.  The  smallest  gain  and  phase  margins 
for  the  controller  are  1 1  dB  and  5  degrees.  Although  the  open  loop  system  is  stable  at 
this  spin  speed  (250  Hz),  these  margins  will  decrease  even  further  as  the  rotor  spin  speed 
increases,  leading  to  instability. 


-600  -500  -400  -300  -200  -100  0 

Phase  [deg] 

Figure  56:  Open  Loop  Gain  vs.  Phase  at  250  Hz, 

Internal  Viscous  Damping  r|v  =  2xl0'6 

Figure  57  presents  the  open  and  closed  loop  frequency  response  of  the  controlled 
bearing  furthest  away  from  the  motor.  The  dashed  line  corresponds  to  the  uncontrolled 
rotor  (open  loop  plant  and  controller  response  combined  from  Figures  55  and  56)  and  the 
solid  line  shows  the  controlled  rotor.  The  difference  between  the  two  represents  the 
action  of  the  magnetic  bearings  upon  the  rotor.  Again,  note  the  sudden  phase  angle 
changes. 
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Plant  and  Controller  Combined  Open/Closed-Loop  FRF:  [PID  kp=. 05  ki=3.00  kd=0.002  pdhz=1000]  Gain  =  1 


Frequency  [Hz] 


Figure  57:  Open  and  Closed  Loop  Frequency  Response  at  250  Hz, 
Internal  Viscous  Damping  t|v  =  2xl0'6 


Based  on  the  results  from  Figures  54-57,  and  based  on  the  earlier  results  of  the 
root  loci  plots  (Figure  51)  it  appears  that  a  PID  controller  is  not  sufficient  to  control  a 
rotor  with  internal  viscous  damping  above  its  critical  speed. 

If  only  internal  hysteretic  damping  was  included  in  the  shaft,  the  simulation 
remained  stable  for  all  spin  speeds  tested.  Figure  58  shows  the  simulated  response  of  the 
rotor  during  a  0  to  400  Hz  ramp  in  one  second.  This  figure  agrees  with  the  FEM  results 
(Table  7)  and  also  is  consistent  with  the  characteristics  of  a  system  with  internal 
hysteretic  damping.  In  this  case,  it  appears  that  the  external  damping  provided  by  the 
magnetic  bearings  was  sufficient  to  overcome  the  amount  of  internal  hysteretic  damping 
in  this  speed  range. 
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simb4.  Rotor  4,  hysteretic  internal  damping,  0-400  Hz  Ramp.  25-Mar-99 


Bearing  1  -  mils  Bearing  2  -  mils 

Figure  58:  Simulation  Prediction  of  Rotor  Response  during  a  0  -  400  Hz  Ramp 
Internal  Hysteretic  Damping  T|h  =  2xl0'4 

Figure  59  shows  the  frequency  response  of  the  rotor  located  at  bearing  number 
two  in  terms  of  magnitude  and  phase  for  a  rotor  speed  of  250  Hz.  The  first  response 
shown  is  the  plant.  The  second  plots  show  the  controller  response  while  the  final  plots 
shows  the  plant  and  controller  combined  response. 
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Figure  59:  (A)  Plant  Response,  (B)  Controller  Response,  (C)  Plant  and 


Controller  Combined  (Open-Loop)  Response  for  a  Rotor  Spin  Speed  of  250  Hz 

Internal  Hysteretic  Damping  T|h  =  2xl0'4 
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Figure  60  plots  the  magnitude  versus  phase  from  the  open  loop  plant  and 
controller  combined  case  shown  in  Figure  59.  This  plot  indicates  that  the  combined  open 
loop  system  is  unstable  at  this  speed,  and  is  consistent  with  the  root  loci  plots  in  Figure 
52. 
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Figure  60:  Open  Loop  Gain  vs.  Phase  at  250  Hz, 

Internal  Hysteretic  Damping  T|h  =  2xl0‘4 

Figure  61  presents  the  open  and  closed  loop  frequency  response  of  the  controlled 
bearing  furthest  away  from  the  motor.  The  dashed  line  corresponds  to  the  uncontrolled 
rotor  (open  loop  plant  and  controller  response  combined  from  Figures  59  and  60)  and  the 
solid  line  shows  the  controlled  rotor. 
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Plant  and  Controller  Combined  Open/Closed-Loop  FRF:  [PID  kp=. 05  ki=3.00kd=0.002pdhz=1000]  Gain  =  1 


Figure  61:  Open  and  Closed  Loop  Frequency  Response  at  250  Hz, 
Internal  Hysteretic  Damping  t|h  =  2x1  O'4 


Although  the  results  from  Figures  58-61  indicate  that  a  PID  controller  can  maintain 
stability  for  the  speed  range  tested,  in  reality  the  inherent  limitations  of  a  PID  controller 
will  eventually  cause  the  rotor  to  become  unstable.  Eventually  as  the  speed  increases,  the 
PID  controller  will  have  to  lower  its  gain  (including  lowering  the  amount  of  external 
damping  applied  to  the  system).  During  this  “roll  off’  period,  the  amount  of  external 
damping  will  fall  below  the  amount  of  internal  hysteretic  damping  and  will  cause  the 
rotor  to  become  unstable,  [21]. 

When  internal  viscous  and  hysteretic  damping  was  added  to  the  shaft,  the  simulation 
predicted  that  the  rotor  would  be  stable  below  its  closed  loop  critical  speed  of  289  Hz. 
Similar  to  the  internal  viscous  damping  case,  this  instability  threshold  is  well  below  the 
FEM  predicted  speed  of  422  Hz  found  in  Table  7.  Figure  62  shows  the  simulated 
response  of  the  rotor  during  a  0  to  400  Hz  ramp  in  one  second. 
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Figure  62:  Simulation  Prediction  of  Rotor  Response  during  a  0  -  400  Hz  Ramp, 
Internal  Viscous  Damping  r|v  =  2x1  O'6,  Internal  Hysteretic  Damping  r|H  =  2xl0'4 


Figure  63  shows  the  frequency  response  of  the  rotor  located  at  bearing  number 
two  in  terms  of  magnitude  and  phase  for  a  rotor  speed  of  250  Hz.  The  first  response 
shown  is  the  plant.  The  second  pair  of  plots  show  the  controller  response  while  the  final 
pair  show  the  plant  and  controller  combined  response.  Note  the  similarities  in  the  phase 
angle  to  Figure  55  (internal  viscous  damping  only). 
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Figure  63:  (A)  Plant  Response,  (B)  Controller  Response,  (C)  Plant  and  Controller 
Combined  (Open-Loop)  Response  for  a  Rotor  Spin  Speed  of  250  Hz, 
Internal  Hysteretic  Damping  T)h  =  2xl0‘4,  Internal  Viscous  Damping  T|v  =  2xl0'6 


Again,  a  number  of  different  PID  controllers  were  tested  in  an  attempt  to  stabilize  the 
rotor  above  its  critical  speed.  Increasing  (or  decreasing)  the  proportional,  integral,  and 
derivative  constants  had  little  effect  on  the  stability  of  the  system.  The  Nichols  plot  in 
Figure  64  shows  how  the  PID  controller  was  limited  by  the  complex  rotordynamics.  The 
gain  and  phase  margins  for  the  controller  are  1 1  dB  and  5  degrees.  Although  the  open 
loop  system  is  stable  at  this  spin  speed  (250  Hz),  these  margins  will  decrease  even  further 
as  the  rotor  spin  speed  increases,  leading  to  instability. 
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Figure  64:  Open  Loop  Gain  vs.  Phase  at  250  Hz,  Internal  Hysteretic  Damping 
T|h  =  2xl0‘4,  Internal  Viscous  Damping  r|v  =  2xl0'6 


Figure  65  presents  the  open  and  closed  loop  frequency  response  of  the  controlled 
bearing  furthest  away  from  the  motor.  The  dashed  line  corresponds  to  the  uncontrolled 
rotor  (open  loop  plant  and  controller  response  combined  case  from  Figures  63  and  64) 
and  the  solid  line  shows  the  controlled  rotor. 


Plant  and  Controller  Combined  Open/Clos ed-Loop  FRF:  [PID  kp=  05  ki=3.00  kd=0.002  pdhz=1000J  Gain  =  1 


Figure  65:  Open  and  Closed  Loop  Frequency  Response,  Internal  Hysteretic 
Damping  r|H  =  2xl0'4,  Internal  Viscous  Damping  r|v  =  2xl0'6 
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Based  on  the  results  from  Figures  63-65,  and  based  on  the  earlier  results  of  the 
root  loci  plots  (Figure  51)  it  appears  that  a  PID  controller  is  not  sufficient  to  control  a 
rotor  with  internal  viscous  and  hysteretic  damping  above  its  critical  speed. 

Again,  it  is  important  to  point  out  how  similar  these  results  are  to  the  case  with 
only  internal  viscous  damping.  This  tends  to  suggest  that  the  stabilizing  effects  of 
internal  viscous  damping  at  speeds  below  its  critical  speed  was  sufficient  to  overcome  the 
destabilizing  effects  of  the  internal  hysteretic  damping  in  the  rotor.  This  is  consistent 
with  the  root  loci  plots  in  Figure  53,  the  stabilizing  term  in  equation  (72),  and  is  also 
shown  in  Figure  16. 

To  summarize,  several  PID  controllers  were  designed  and  implemented  in  an 
effort  to  stabilize  a  rotor  system  with  internal  damping.  These  efforts  focused  on 
controller  design  in  the  frequency  domain  with  subsequent  testing  in  the  full  simulation. 
For  the  cases  with  internal  viscous  damping,  the  PID  controller  was  able  to  maintain 
stability  for  sub-critical  speeds  only.  For  the  internal  hysteretic  damping  case,  the  PID 
controller  was  able  to  maintain  control  for  all  speeds  tested,  although  in  reality  the 
inherent  limitations  of  a  PID  controller  will  eventually  cause  the  rotor  to  become 
unstable.  Table  8  compares  the  simulation  and  FEM  results  for  all  cases  of  internal 
damping. 

Although  PID  controllers  were  effective  in  stabilizing  the  rotor  systems  at  sub- 
critical  speeds,  the  model  developed  during  this  thesis  showed  that  these  controllers  were 
unable  to  counteract  the  destabilizing  effects  of  internal  damping  during  supercritical 
operation.  Other  more  complex  controllers  should  be  designed  specifically  to  achieve 
improved  control  and  increase  the  instability  threshold  for  the  system. 
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Other  important  findings  regarding  internal  damping  include: 

-  The  erratic  changes  in  the  phase  angle  for  systems  with  internal  damping  could  be 
a  factor  in  determining  the  stability  of  the  system.  Further  study  should  be 
performed  to  determine  specifically  how  the  phase  angle  and  internal  damping  are 
related. 

-  The  results  of  the  rotor  with  internal  viscous  damping  were  very  similar  to  the 
results  of  the  case  where  both  internal  viscous  and  hysteretic  damping  were 
included.  This  tends  to  suggest  that  the  stabilizing  effects  of  internal  viscous 
damping  at  speeds  below  its  critical  speed  was  sufficient  to  overcome  the 
destabilizing  effects  of  the  internal  hysteretic  damping  in  the  rotor. 

Table  8:  Comparison  of  Instability  Threshold  Predictions 
FEM  Model:  Bearing  Stiffness  =  50  lbs/in,  Bearing  Damping  =  1.5xl0'2  lbs-sec/in 
Simulation:  Kp  =  0.05  lb/mil,  Ki  =  3.00  lb/mil-s,  Kd  =  0.002  lb-sec/mil, 
Derivative  Bandwidth  =  1000  Hz 


Internal  Damping 

Critical  Speeds  (Hz) 

Instability  T1 

hreshold  (Hz) 

■m 

IBH 

T|V  =  2x1 0'6 

272.5 

289 

441 

289 

T|h  =  0.0002 

272.5 

289 

Stable  for 
all  Speeds 

Stable  for  all 
Speeds 

T|v  =  2x10  T|h  =  0.0002 

272.5 

289 

425 

289 

1  Extended  FEM  model  has  constant  bearing  properties 

2  Modified  Simulation 
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Chapter  8  Summary  and  Conclusions 


In  the  first  few  chapters  of  this  thesis  the  concepts  of  external  and  internal 
damping  were  explored.  Internal  damping  was  defined  and  divided  into  two  distinct 
types:  internal  viscous  and  internal  hysteretic  damping.  Next,  internal  damping  was 
included  in  a  rotating  system  and  it  was  shown  how  under  certain  conditions,  internal 
damping  can  cause  the  rotor  to  become  unstable.  The  following  general  conclusions 
regarding  internal  and  external  damping  and  their  effects  on  stability  were  made: 

-  If  the  external  and  internal  damping  are  both  viscous  in  nature,  then  the  rotor  will 
always  be  stable  in  its  sub-critical  speed  range.  Stability  can  be  extended  into  the 
supercritical  regime  by  adding  external  damping  according  to  the  results  obtained  by 
Smith,  equation  (56). 

-  If  the  external  and  internal  damping  are  both  hysteretic,  the  rotor  will  be  unstable  for 
all  speeds  unless  the  external  damping  is  large  enough  to  establish  stability.  In  this 
case,  the  rotor  will  be  stable  for  all  speeds. 

-  If  the  external  damping  is  viscous  and  the  internal  damping  is  hysteretic,  the  rotor 
will  be  unstable  in  a  range  of  speeds  beginning  at  £2  =  0+.  As  the  external  viscous 
damping  increases,  this  speed  range  shrinks  toward  £2  =  0  until  the  external  damping 
reaches  a  value  that  results  in  stability  for  all  speeds,  i.e.  for  Q  >  0. 

After  the  general  characteristics  of  internal  damping  were  discussed,  a  finite 
element  model  that  included  internal  damping  was  developed  and  validated  using  a 
standard  test  case  from  the  references.  This  model  was  then  used  to  produce  state  space 
matrices  that  fully  described  a  rotor  with  and  without  internal  damping.  To  make  use  of 
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these  matrices,  a  Matlab™  magnetic  bearing  simulation  developed  by  Draper  Laboratory 
(Scholten)  was  modified  to  include  internal  damping.  The  modified  simulation  was  able 
to  predicted  rigid  body  modes,  critical  speeds,  rotor  displacements  and  instability 
thresholds. 

Prior  to  entering  the  rotor  model  with  internal  damping  into  the  simulation,  the 
original  FEM  model  (Antkowiak)  and  simulation  (Scholten)  predictions  were  compared 
to  actual  rotor  test  results.  The  FEM  model  critical  speed  prediction  was  slightly  lower 
than  the  actual  rotor  rig  critical  speed.  This  was  probably  due  to  the  simplicity  of  the 
FEM  model  since  it  did  not  take  into  account  rotor  unbalance  and  it  did  not  model  any 
other  overall  system  effects  such  as  magnetic  bearing  material  non-linearity,  actuator  lag, 
etc. 

In  contrast  to  the  simple  FEM  model,  the  high  fidelity  simulation  predicted  a  first 
critical  speed  that  closely  matched  the  actual  rotor  test  data,  although  the  rigid  body  mode 
predictions  were  not  as  accurate.  With  regard  to  the  displacements  predicted  by  the 
simulation,  good  agreement  was  seen  in  terms  of  the  predicted  and  actual  displacements 
for  the  bearing  furthest  away  from  the  elastic  coupling,  although  the  simulation  predicted 
somewhat  lower  displacements  for  the  closer  bearing.  This  suggests  that  more  work 
needs  to  be  done  on  the  elastic  coupling  block  of  the  simulation.  In  terms  of  the  vertical 
displacements  at  the  bearings,  the  simulation  had  good  correlation  to  the  actual  rotor  data, 
while  the  predicted  horizontal  displacements  were  slightly  smaller  than  the  actual  rotor 
data. 

The  only  area  where  the  simulation  had  poor  accuracy  with  regard  to  the  actual 
rotor  displacements  was  in  the  speed  range  near  the  rigid  body  mode.  Since  the 
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simulation  predicted  the  rigid  body  mode  at  47  Hz,  its  displacement  predictions  at  the 
actual  rigid  body  mode  of  52  Hz  would  not  be  as  large. 

After  verifying  the  high  fidelity  simulation  and  incorporating  the  destabilizing 
effects  of  internal  damping,  several  PID  controllers  were  designed  and  implemented  in  an 
effort  to  stabilize  a  rotor  system  with  internal  damping.  These  efforts  focused  on 
controller  design  in  the  frequency  domain  with  subsequent  testing  in  the  full  simulation. 
For  the  cases  with  internal  viscous  damping,  the  PID  controller  was  able  to  maintain 
stability  for  sub-critical  speeds  only.  For  the  internal  hysteretic  damping  case,  the  PID 
controller  was  able  to  maintain  control  for  all  speeds  tested,  although  in  reality  the 
inherent  limitations  of  a  PID  controller  will  eventually  cause  the  rotor  to  become  unstable 
above  its  critical  speed.  Other  important  findings  regarding  internal  damping  and  the  PID 
controller  included: 

-  The  erratic  changes  in  the  phase  angle  for  systems  with  internal  damping  could  be 
a  factor  in  determining  the  stability  of  the  system.  Further  study  should  be 
performed  to  determine  specifically  how  the  phase  angle  and  internal  damping  are 
related. 

-  The  results  of  the  rotor  with  internal  viscous  damping  were  very  similar  to  the 
results  of  the  case  where  both  internal  viscous  and  hysteretic  damping  were 
included.  This  tends  to  suggest  that  the  stabilizing  effects  of  internal  viscous 
damping  at  speeds  below  its  critical  speed  was  sufficient  to  overcome  the 
destabilizing  effects  of  the  internal  hysteretic  damping  in  the  rotor. 
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To  continue  this  study,  a  number  of  tasks  should  be  undertaken.  First,  a  proper 
measurement  of  the  internal  viscous  and  hysteretic  damping  of  the  shaft  should  be  made. 
Next,  a  higher  speed  motor  must  be  acquired  to  test  where  the  rotor  actually  becomes 
unstable.  Based  on  these  actual  rotor  results,  the  internal  damping  models  might  need  to 
be  extended.  In  addition,  instability  threshold  predictions  of  shaft  with  an  additional 
inertia,  such  as  a  disc,  should  be  compared  to  actual  test  data.  This  would  lead  to  a  more 
useful  simulation,  since  it  is  more  representative  of  real-world  applications.  A  final 
suggestion  for  further  study  in  the  realm  of  magnetic  bearings  and  controls  should 
include  testing  whether  centralized  control,  as  opposed  to  the  current  decentralized 
control,  could  improve  the  system  stability  margin. 

The  PID  controllers  were  effective  in  stabilizing  the  rotor  systems  at  sub-critical 
speeds.  However,  the  model  developed  during  this  thesis  showed  that  these  controllers 
were  unable  to  counteract  the  destabilizing  effects  of  internal  damping  during 
supercritical  operation.  Other  more  complex  controllers  should  be  designed  and  tested 
specifically  to  achieve  improved  control  and  increase  the  instability  threshold  for  the 
system.  This  improved  rotordynamic  model  and  magnetic  bearing  simulation  is  now 
available  to  test  these  designs  in  the  supercritical  regime. 
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Appendix  A  -  Detailed  Information  on  Draper  Magnetic  Bearing  Test  Apparatus 


Base 

Figure  Al:  Draper  Magnetic  Bearing  Test  Rig 

Motor  Base  - 

The  base  structure  for  the  test  rig  is  from  a  Bently  Nevada  RK-3  rotor  kit. 
This  kit  includes  an  electric  motor  mounted  on  a  base  plate  that  drives  the  shaft.  The 
motor  and  controller  have  the  capability  to  operate  in  the  range  of  0  to  10,000  rpm.  The 
motor  is  connected  to  the  shaft  via  an  elastic  coupling. 

Rotor  - 

The  rotating  assembly  is  comprised  of  two  actuator  rotors  designed  by  Draper 
Laboratory  pressed  onto  a  3/8”  diameter  stainless  steel  shaft.  The  actuator  is  an 
aluminum  spool  surrounded  by  multiple  magnetic  material  laminations  and  a  ferrous  iron 
ring,  which  is  used  as  the  sensor  plate  for  the  eddy  current  sensors.  A  threaded  brass  nut 
is  used  to  secure  the  laminations  to  the  spool.  The  iron  rig  is  mounted  with  a  shrink  fit. 
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Magnetic  Actuators- 

A  four-pole  hetro  polar  bearing  design  at  Draper  Laboratory  was  used  on  the  rotor 
test  rig.  A  frameless  stator  design  was  selected  so  that  the  stator  laminations  could  easily 
be  mounted  in  a  common  aluminum  housing.  A  nominal  radial  air  gap  of  0.012”  was 
selected. 

The  lamination  stack  for  the  stator  and  rotor  were  made  out  of  Carpenter  Hiperco 
alloy  50  in  14  mil  sheets.  Hiperco50  is  an  iron-cobalt-vandium  soft  magnetic  alloy.  This 
alloy  was  advertised  to  have  a  saturation  of  2  Tesla  when  heat-treated  accordingly. 

The  stator  was  designed  to  accommodate  120  turns  of  #20  gauge  wire.  To 
minimize  thermal  heating  in  the  copper  coil,  a  maximum  current  density  of  5,000 
amps/in2  was  used.  For  desirable  operation  in  the  linear  range  of  the  BH  curve,  less  that 
1.74  tesla,  and  a  bearing  load  capacity  of  10  lbs,  sized  the  cross  sectional  area  of  the  pole 
to  be  0.025  in2.  This  calculation  dictated  that  a  stack  of  9  laminations  (0.126  in  thick  by 
0.2  pole  width)  was  used  in  the  stator  design,  [21]. 

Touchdown  Bearings- 

In  case  of  excess  shaft  vibration  or  controller  failure,  touchdown  bearings  were 
provided  to  protect  the  stators  of  the  magnetic  bearings.  They  provided  a  radial  clearance 
of  0.008”  from  the  center  of  the  bearing. 


A-2 


Computer  hardware  and  software- 

The  digital  signal  processing  (DSP)  hardware  is  based  on  a  single  Texas  Instruments 
C40  DSP  operating  at  50  megahertz.  Several  PC-hosted  circuit  boards  were  used  for 
processing,  digital  to  analog  conversion,  and  analog  to  digital  conversions.  All  of  these 
boards  were  manufactured  by  dSpace,  inc.  Along  with  their  hardware,  dSpace  provided 
low-level  support  software  libraries  for  running  and  monitoring  the  DSP.  The  magnetic 
bearing  control  software  running  on  the  DSP  was  written  by  Draper  Laboratory  in  the  C 
language  (Scholten).  For  speed  and  calculation  robustness,  linear  transfer  functions  were 
coded  in  assembly  language  as  Z-transform  biquad  filters,  (Scholten,  [22]).  The  sample 
rate  was  10  kHz. 

The  DSP  system  was  hosted  on  a  Pentium  PC.  It  provided  support  for  the 
following: 

-  Data  logging,  using  Trace40  application  from  dSpace; 

-  Analysis  of  logged  data,  using  Draper  developed  software  routines; 

-  System  Simulation  using  the  Matlab  and  Simulink  Applications; 

-  A  custom  user  interface  written  in  C  which  provided  for  monitoring  the 
dSpace  real  time  (Scholten). 
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Appendix  B  -  A  Timoshenko  Rotor  Element 
No  Internal  Damping 

In  line  with  its  work  on  magnetic  bearings,  The  Charles  Stark  Draper  Laboratory 
has  produced  a  rotor  dynamic  finite  element  model  that  is  based  on  Timoshenko  beam 
theory.  A  Timoshenko  rotor  uses  an  element  with  2  nodes  and  8  Degrees  of  Freedom  and 
uses  third  order  shape  functions  to  describe  the  bending  of  the  element.  All  additional 
inertia  are  assumed  to  be  rigid  discs  with  lump  mass  properties  and  the  bearings  are 
assumed  to  be  discrete  and  linear.  The  model  includes  rotary  inertia,  gyroscopic 

moments,  and  shear  deformation  effects  (<I> ).  Following  the  derivation  described  by 
Nelson  [16, 17]  the  shape  functions  are: 


Translational  Shape  Functions: 


Rotational  Shape  Functions: 


(B.la,b,c,d) 
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0,  =~-rk+®v,] 

1  +  <I> 

„  1  [if  Js\  Is 

0  —  — —  —  —  6  —  +  6  — 

l  +  0>  /  /  / 


©  ,=-L  1-41-1+3141  +  ®(l-4 


/  / 


©,=-!-  i  6141-614 


(B.2a,b,c,d) 


l  +  0>  l  Z  J  / 


e,  =-X=  -2|- 1+3141  +®i 


i+®  /  / 


The  matrix  [A]  transforms  the  nodal  displacement  vector  (V,  W)  into  the  rotor 
translational  displacements  (Vbend,  Wbend),  while  the  matrix  [©]  converts  the 
displacements  into  the  rotor  cross-section  center  line  translations  (VShear,  WShear)- 


r  1  A,  0  0  A2  A3  0  0  A4 

U=[o  A,  -A2  0  0  A3  -A4  0 

r  l_r  0  -©!  02  0  0  -03  ©4  0 

0,  0  0  ©2  ©3  0  0  04 


(B.3a,b) 


From  the  rotor  configuration  shown  in  Figure  Bl,  the  rotations  B  and  T  are  defined  as: 
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Figure  Bl:  Rotor  System  Configuration,  Nelson 

Using  the  Lagrangian  approach,  Nelson  calculated  the  potential  and  kinetic 
energy  of  the  rotating  element.  The  rotor  potential  energy  of  the  rotor  was  stored  in  two 
forms,  bending  deformation,  and  shear  deformation  (axial  loads  neglected): 


(B.6) 


where  the  shear  deformation  factor  is  given  by  K  and  the  second  derivative  of  the  bending 
deformation  in  the  y  direction  is  represented  by  V*end .  The  rotor  kinetic  energy  included 
Timoshenko  effects  of  rotary  inertia,  shear  deformation,  and  the  gyroscopic  energy: 


r  1 


(B.7) 


+  J-  pIpQ2ds  -  J  IpQj3Tdz 


Nelson,  then  proceeded  to  derive  the  following  matrices: 
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m2 

0  m2 


where  the  elements  in  the  translational  mass  matrix  are  defined  by: 


m,  =156  +  2490 +  140O2 
m2=l2( 4+70  +  3.502) 
m3  =/(22  + 38.50 +17.502) 
m4  =  54 +  1260  +  700 2 
m5  =/(13  +  31.50  +  17.502) 
m6=l2( 3  +  70  +  3.502) 


(B.9a,b,c,d,e,f) 
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where  the  elements  in  the  rotational  mass  matrix  are  defined  by: 
n,  =36 

n2  =/2(4  +  5O  +  10O2) 
n3  =/(3-150) 
n4=/2(l  +  50-502) 


(B.lla,b,c,d) 
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where  the  individual  elements  in  the  stiffness  matrix  are  defined  by: 

*,=12 

k2  =  ( length  l)2 (4  +  O) 
k3  =  61 

k4=l2(  2-0) 


(B.13,a,b,c,d) 
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where  the  individual  element  in  the  gyroscopic  matrix  are  defined  by: 


S,=36 

g2  =(4  +  5^+10O2)/ 
g3=  (-3  +  15*)/ 
§4=(-l-50  +  502)/2 


(B.15a,b,c,d) 


By  the  application  of  Hamilton’s  extended  principle,  and  using  the  above  matrices 
and  energy  and  work  functions,  Nelson  produced  the  following  undamped  matrix 
equations  of  motion: 
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m,raJHMrJ){q}-Q[G]{q}  +  [Kbend]{q}  =  {f} 


(B.16) 


where: 

{q}  -  fixed  frame  physical  coordinates 
{f}  -  fixed  frame  external  forces 
Q  -  rotor  spin  speed 

The  eigenvalues  for  the  undamped  equation  of  motion  occur  in  conjugate  pairs 
and  are  given  by:  Sj  =  ±iCQj ,  where  the  imaginary  pair  ±co  represents  the  forward  and 

backward  whirl  frequencies  of  the  shaft. 


Addition  of  Internal  Damping 


Based  on  the  work  of  Zorzi  and  Nelson,  [2],  the  previous  derivation  of  the 
rotordynamic  equations  of  motion  can  be  modified  to  include  the  contributions  of  internal 
damping.  Nelson  used  both  of  the  internal  viscous  and  hysteretic  damping  models 
discussed  in  the  previous  chapters.  Adding  this  combination  of  internal  damping  terms 
into  the  constitutive  relationship  yields: 

f  f  ^  j 


a  =E< 


Vv  + 


W1  +  77« 


(B.17) 


Further,  by  defining  the  whirl  radius  as  R  and  the  shaft  radius  as  r  (see  Figure  Bl),  the 
strain  and  strain  rate  are: 
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(B.18a,b) 


ot  o2x 


Again  from  Figure  Bl,  the  internal  bending  moments  Mz  and  My  can  be  expressed  as: 

My  =\[w  +  rsin(£lt)]jxdr(rd(Q.t)) 

*  _  .  (B.19atb) 

Mz  =  j  -  [V  +  rcos(£lt)}jxdr(rd(Q.t)) 

A 


After  substituting  the  modified  constitutive  relationship  into  the  internal  bending 
moments  and  performing  the  integrations,  the  moments  can  be  expressed  as: 
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Now  placing  the  energy  contributions  from  these  moment  equations  into  the 
appropriate  kinetic  or  potential  energy  equations  defined  in  section  4.1,  gives: 


(B.22) 
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Solving  the  equations  of  motion  with  internal  damping  results  in  eigenvalues  in 
the  form  of:  Sj  =  Xs  ±iC0j ,  where  ©  provides  shaft  whirl  and  X  provides  an  orbit 

growth/decay  after  a  perturbation. 
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Appendix  C:  Input  Files  for  Draper  Rotor  FEM  Program 
Input  File  for  Rotor  -  No  Internal  Damping,  No  Constraints 


ELEMENT 

1 

3  0.500 

0.3125 

ELEMENT 

2 

1  0.675 

0.3748 

ELEMENT 

3 

1  0.675 

0.3748 

ELEMENT 

4 

2  0.360 

0.3748 

ELEMENT 

5 

2  0.500 

0.3748 

ELEMENT 

6 

2  0.440 

0.3748 

ELEMENT 

7 

1  1.000 

0.3748 

ELEMENT 

8 

1  1.000 

0.3748 

ELEMENT 

9 

1  1.000 

0.3748 

ELEMENT 

10 

1  1.000 

0.3748 

ELEMENT 

11 

1  1.000 

0.3748 

ELEMENT 

12 

1  1.000 

0.3748 

ELEMENT 

13 

1  1.000 

0.3748 

ELEMENT 

14 

1  1.000 

0.3748 

ELEMENT 

15 

1  1.000 

0.3748 

ELEMENT 

16 

2  0.440 

0.3748 

ELEMENT 

17 

2  0.500 

0.3748 

ELEMENT 

18 

2  0.360 

0.3748 

ELEMENT 

19 

1  0.625 

0.3748 

ELEMENT 

20 

1  0.625 

0.3748 

MATERIAL  1 

29.5E+06  .33  7.27-04 

MATERIAL  2 

69.3E+06  .33  8.40-04 

MATERIAL  3 

10.0E+06  .45  7.50-04 

CONSTRAINT  1 

000.0  000.0 

CONSTRAINT  5 

0.000E+00  0.000E+00 

CONSTRAINT  6 

0.000E+04  0.000E+04 

CONSTRAINT  8 

0.000E+00  0.000E+00 

CONSTRAINT  11 

0.O00E+00  0.000E+00 

CONSTRAINT  15 

O.OOOE+OO  0.000E+00 

CONSTRAINT  17 

0.000E+04  0.000E+04 

CONSTRAINT  18 

0.000E+00  0.000E+00 

CONSTRAINT  21 

O.OOOE+OO  0.000E+00 

ADDMASS  4 

32.00E-06  1.500E-06  2.700E-06 

ADDMASS  5 

68.00E-06  4.100E-06  5.800E-06 

ADDMASS  6 

62.00E-06  3.600E-06  5.300E-06 

ADDMASS  7 

26.00E-06  1.200E-06  2.200E-06 

ADDMASS  16 

26.00E-06  1.200E-06  2.200E-06 

ADDMASS  17 

62.00E-06  3.600E-06  5.300E-06 

ADDMASS  18 

68.00E-06  4.100E-06  5.800E-06 

ADDMASS  19 

32.00E-06  1.500E-06  2.700E-06 
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Input  File  for  Rotor  -  No  Internal  Damping,  With  Constraints  from  Estimated 
Bearing 


ELEMENT 

1 

3  0.500 

0.3125 

ELEMENT 

2 

1  0.675 

0.3748 

ELEMENT 

3 

1  0.675 

0.3748 

ELEMENT 

4 

2  0.360 

0.3748 

ELEMENT 

5 

2  0.500 

0.3748 

ELEMENT 

6 

2  0.440 

0.3748 

ELEMENT 

7 

1  1.000 

0.3748 

ELEMENT 

8 

1  1.000 

0.3748 

ELEMENT 

9 

1  1.000 

0.3748 

ELEMENT 

10 

1  1.000 

0.3748 

ELEMENT 

11 

1  1.000 

0.3748 

ELEMENT 

12 

1  1.000 

0.3748 

ELEMENT 

13 

1  1.000 

0.3748 

ELEMENT 

14 

1  1.000 

0.3748 

ELEMENT 

15 

1  1.000 

0.3748 

ELEMENT 

16 

2  0.440 

0.3748 

ELEMENT 

17 

2  0.500 

0.3748 

ELEMENT 

18 

2  0.360 

0.3748 

ELEMENT 

19 

1  0.625 

0.3748 

ELEMENT 

20 

1  0.625 

0.3748 

MATERIAL 

1 

29.5E+06  .33 

7.27-04 

MATERIAL 

2 

69.3E+06  .33 

8.40-04 

MATERIAL 

3 

10.0E+06  .45 

7.50-04 

CONSTRAINT 

1  000.0 

000.0 

CONSTRAINT 

5  0.000E+00  0.000E+00 

CONSTRAINT  6  5.000E+01  5.000E+01  1.500E-02  1.500E-02 

CONSTRAINT  8  O.OOOE+OO  0.000E+00 

CONSTRAINT  11  0.000E+00  0.000E+00 

CONSTRAINT  15  0.000E+00  0.000E+00 

CONSTRAINT  17  5.000E+01  5.000E+01  1.500E-02  1.500E-02 

CONSTRAINT  18  0.000E+00  0.000E+00 

CONSTRAINT  21  0.000E+00  0.000E+00 

ADDMASS  4  32.00E-06 1.500E-06  2.700E-06 

ADDMASS  5  68.00E-064.100E-06  5.800E-06 

ADDMASS  6  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  7  26.00E-06 1.200E-06  2.200E-06 

ADDMASS  16  26.00E-06  1.200E-06  2.200E-06 

ADDMASS  17  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  1 8  68.00E-06  4. 1 00E-06  5.800E-06 

ADDMASS  19  32.00E-06  1 .500E-06  2.700E-06 
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Input  File  for  Rotor  -  Internal  Damping  Included,  No  Constraints 


ELEMENT 

1 

3  0.500 

0.3125 

ELEMENT 

2 

1  0.675 

0.3748 

ELEMENT 

3 

1  0.675 

0.3748 

ELEMENT 

4 

2  0.360 

0.3748 

ELEMENT 

5 

2  0.500 

0.3748 

ELEMENT 

6 

2  0.440 

0.3748 

ELEMENT 

7 

1  1.000 

0.3748 

ELEMENT 

8 

1  1.000 

0.3748 

ELEMENT 

9 

1  1.000 

0.3748 

ELEMENT 

10 

1  1.000 

0.3748 

ELEMENT 

11 

1  1.000 

0.3748 

ELEMENT 

12 

1  1.000 

0.3748 

ELEMENT 

13 

1  1.000 

0.3748 

ELEMENT 

14 

1  1.000 

0.3748 

ELEMENT 

15 

1  1.000 

0.3748 

ELEMENT 

16 

2  0.440 

0.3748 

ELEMENT 

17 

2  0.500 

0.3748 

ELEMENT 

18 

2  0.360 

0.3748 

ELEMENT 

19 

1  0.625 

0.3748 

ELEMENT 

20 

1  0.625 

0.3748 

MATERIAL  1  29.5E+06  .33  7.27-04  2.00E-04  2.00E-06 

MATERIAL  2  69.3E+06  .33  8.40-04  2.00E-04  2.00E-06 

MATERIAL  3  10.0E+06  .45  7.50-04 

CONSTRAINT  1  000.0  000.0 

CONSTRAINT  5  0.000E+00  0.000E+00 

CONSTRAINT  6  0.000E+04  0.000E+04 
CONSTRAINT  8  0.000E+00  0.000E+00 

CONSTRAINT  11  0.000E+00  0.000E+00 

CONSTRAINT  15  0.000E+00  0.000E+00 

CONSTRAINT  17  0.000E+04  0.000E+04 

CONSTRAINT  18  0.000E+00  0.000E+00 

CONSTRAINT  21  0.000E+00  0.000E+00 

ADDMASS  4  32.00E-06  1 .500E-06  2.700E-06 

ADDMASS  5  68.00E-06  4.100E-06  5.800E-06 

ADDMASS  6  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  7  26.00E-06  1 .200E-06  2.200E-06 

ADDMASS  16  26.00E-06  1.200E-06  2.200E-06 

ADDMASS  17  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  18  68.00E-06  4.100E-06  5.800E-06 

ADDMASS  19  32.00E-06  1.500E-06  2.700E-06 
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Input  File  for  Rotor  -  Internal  Damping  Included,  With  Constraints 


ELEMENT 

1 

3  0.500 

0.3125 

ELEMENT 

2 

1  0.675 

0.3748 

ELEMENT 

3 

1  0.675 

0.3748 

ELEMENT 

4 

2  0.360 

0.3748 

ELEMENT 

5 

2  0.500 

0.3748 

ELEMENT 

6 

2  0.440 

0.3748 

ELEMENT 

7 

1  1.000 

0.3748 

ELEMENT 

8 

1  1.000 

0.3748 

ELEMENT 

9 

1  1.000 

0.3748 

ELEMENT 

10 

1  1.000 

0.3748 

ELEMENT 

11 

1  1.000 

0.3748 

ELEMENT 

12 

1  1.000 

0.3748 

ELEMENT 

13 

1  1.000 

0.3748 

ELEMENT 

14 

1  1.000 

0.3748 

ELEMENT 

15 

1  1.000 

0.3748 

ELEMENT 

16 

2  0.440 

0.3748 

ELEMENT 

17 

2  0.500 

0.3748 

ELEMENT 

18 

2  0.360 

0.3748 

ELEMENT 

19 

1  0.625 

0.3748 

ELEMENT 

20 

1  0.625 

0.3748 

MATERIAL  1  29.5E+06  .33  7.27-04  2.00E-04  2.00E-06 

MATERIAL  2  69.3E+06  .33  8.40-04  2.00E-04  2.00E-06 

MATERIAL  3  10.0E+06  .45  7.50-04 

CONSTRAINT  1  000.0  000.0 

CONSTRAINT  5  O.OOOE+OO  0.000E+00 

CONSTRAINT  6  5.000E+01  5.000E+01  1.500E-02  1.500E-02 

CONSTRAINT  8  0.000E+00  0.000E+00 

CONSTRAINT  1 1  0.000E+00  0.000E+00 

CONSTRAINT  15  0.000E+00  0.000E+00 

CONSTRAINT  17  5.000E+01  5.000E+01  1.500E-02  1.500E-02 

CONSTRAINT  18  0.000E+00  0.000E+00 

CONSTRAINT  21  0.000E+00  0.000E+00 

ADDMASS  4  32.00E-06  1.500E-06  2.700E-06 

ADDMASS  5  68.00E-06  4.100E-06  5.800E-06 

ADDMASS  6  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  7  26.00E-06  1.200E-06  2.200E-06 

ADDMASS  16  26.00E-06  1 .200E-06  2.200E-06 

ADDMASS  17  62.00E-06  3.600E-06  5.300E-06 

ADDMASS  1 8  68.00E-06  4. 1 00E-06  5.800E-06 

ADDMASS  19  32.00E-06  1.500E-06  2.700E-06 


C-4 


